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Abstract 

This paper develops a general framework for solving a variety of convex cone problems that 
frequently arise in signal processing, machine learning, statistics, and other fields. The approach 
works as follows: first, determine a conic formulation of the problem; second, determine its dual; 
third, apply smoothing; and fourth, solve using an optimal first-order method. A merit of this 
approach is its flexibility: for example, all compressed sensing problems can be solved via this 
approach. These include models with objective functionals such as the total-variation norm, 
||Wa;||i where W is arbitrary, or a combination thereof. In addition, the paper also introduces 
a number of technical contributions such as a novel continuation scheme, a novel approach 
for controlling the step size, and some new results showing that the smooth and unsmoothed 
problems are sometimes formally equivalent. Combined with our framework, these lead to novel, 
stable and computationally efficient algorithms. For instance, our general implementation is 
competitive with state-of-the-art methods for solving intensively studied problems such as the 
LASSO. Further, numerical experiments show that one can solve the Dantzig selector problem, 
for which no efficient large-scale solvers exist, in a few hundred iterations. Finally, the paper is 
accompanied with a software release. This software is not a single, monolithic solver; rather, it is 
a suite of programs and routines designed to serve as building blocks for constructing complete 
algorithms. 

Keywords. Optimal first-order methods, Nesterov's accelerated descent algorithms, proximal 
algorithms, conic duality, smoothing by conjugation, the Dantzig selector, the LASSO, nuclear- 
norm minimization. 



1 Introduction 

1.1 Motivation 

This paper establishes a general framework for constructing optimal first-order methods for solving 
certain types of convex optimization programs that frequently arise in signal and image processing, 
statistics, computer vision, and a variety of other fields^] In particular, we wish to recover an 

1 The meaning of the word 'optimal' shall be made precise later. 
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unknown vector xq £ W 1 from the data y £ M. m and the model 

y = Ax + z; (1.1) 

here, A is a known m x n design matrix and z is a noise term. To fix ideas, suppose we find 
ourselves in the increasingly common situation where there are fewer observations/measurements 
than unknowns, i.e., m < n. While this may seem a priori hopeless, an impressive body of recent 
works has shown that accurate estimation is often possible under reasonable sparsity constraints on 
Xq. One practically and theoretically effective estimator is the Dantzig selector introduced in |15| . 
The idea of this procedure is rather simple: find the estimate which is consistent with the observed 
data and has minimum t\ norm (thus promoting sparsity). Formally, assuming that the columns 
of A are normalized]^] the Dantzig selector is the solution to the convex program 

minimize IMIi 2^ 
subject to \\A*(y- Ax)\\oo < 5, ^ ' ' 

where 5 is a scalar. Clearly, the constraint is a data fitting term since it asks that the correlation 
between the residual vector r = y — Ax and the columns of A is small. Typically, the scalar 5 is 
adjusted so that the true xo is feasible, at least with high probability, when the noise term z is 
stochastic; that is, 5 obeys ||A*£||oo < 5 (with high probability). Another effective method, which 
we refer to as the LASSO [45] (also known as basis pursuit denoising, or BPDN), assumes a different 
fidelity term and is the solution to 

minimize \\ x \\i q > 

subject to \\y — Ax\\2 < e, 

where e is a scalar, which again may be selected so that the true vector is feasible. Both of these 
estimators are generally able to accurately estimate nearly sparse vectors and it is, therefore, of 
interest to develop effective algorithms for each that can deal with problems involving thousands 
or even millions of variables and observations. 



There are of course many techniques, which are perhaps more complicated than (1.2) and (1.3) 



for recovering signals or images from possibly undersampled noisy data. Suppose for instance that 



we have noisy data y (1.1) about an n x n image xq> that is, [xo]ij i s an fi 2 -array of real numbers. 



Then to recover the image, one might want to solve a problem of this kind: 

minimize || Wx||i + A||x||tv (141 
subject to \\y — Ax\\2 < e, 

where W is some (possibly nonorthogonal) transform such as an undecimated wavelet transform 
enforcing sparsity of the image in this domain, and || • ||tv is the total- variation norm defined as 

\\ x \\tv ■= ^2 vVl* + 1 J] ~ x ihj]\ 2 + + 1] - x[i,j]\ 2 . 



The motivation for (1.4) is to look for a sparse object in a transformed domain while reducing 



artifacts due to sparsity constraints alone, such as Gibbs oscillations, by means of the total- variation 



There is a slight modification when the columns do not have the same norm, namely, WD' 1 A* (y - Ax)]^ < 5, 
where D is diagonal and whose diagonal entries are the £2 norms of the columns of A. 
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norm [10 13,43 . The proposal (1.4 1 appears computationally more involved than both (1.2) and 
(1.3), and our goal is to develop effective algorithms for problems of this kind as well. 

To continue our tour, another problem that has recently attracted a lot attention concerns the 
recovery of a low-rank matrix Xq from undersampled data 

y = A(X ) + z, (1.5) 

where A : W niXri2 — > M m is a linear operator supplying information about Xq. An important 
example concerns the situation where only some of the entries of Xq are revealed, A{Xq) = [Ao]ij : 
(i, j) G E C [m] x [712], and the goal is to predict the values of all the missing entries. It has been 



shown [12 16 24 that an effective way of recovering the missing information from y and the model 
(1.5) is via the convex program 

minimize ,~ „s 

subject to X G C. ^ ' 

Here, is the sum of the singular values of the matrix X, a quantity known as the nuclear 

norm of X. (\\X\\* is also the dual of the standard operator norm ||X||, given by the largest singular 
value of X). Above, C is a data fitting set, and might be {X : A(X) = y} in the noiseless case, 
or {X : \\A*(y - ^t(X))||oo < 5} (Dantzig selector-type constraint) or {X : \\y - A(X)\\ 2 < e} 
(LASSO-type constraint) in the noisy setup. We are again interested in computational solutions 
to problems of this type. 

1.2 The literature 

There is of course an immense literature for solving problems of the types described above. Consider 
the LASSO, for example. Most of the works [§[^[22j[27j[39j|48j[5UJ[52] are concerned with the 
unconstrained problem 

minimize g \\Ax — &H2 + A||x||i, (1-7) 
which differs from (|1.3[) in that the hard constraint \\Ax — 6|| 2 < e is replaced with a quadratic 



penalty 1 1| A2; — 6|j|. There are far fewer methods specially adapted to (1.3); let us briefly 



discuss some of them. SPGL1 47 is an excellent solver specifically designed for (1.3). The issue is 



that at the moment, it cannot handle important variations such as 

minimize ||Wx||i 
subject to \\y — Ax\\2 < e, 



where W is a nonorthogonal transform as in (1.4). The main reason is that SPGL1 — as with almost 



all first-order methods for that matter — relies on the fact that the proximity operator associated 
with the t\ norm, 

x(z; t) = argmin |t~ \\x — + ||ar||i, (1-8) 

is efficiently computable via soft-thresholding. This is not the case, however, when ||x||i is replaced 
by a general term of the form ||Wx||i. NESTA [4] can efficiently deal with an objective functional 
of the form [|Wx||i, but it requires repeated projections onto the feasible set; see also [l] for a 
related approach. Hence, NESTA is efficient when AA* is a projector or, more generally, when the 
eigenvalues of AA* are well clustered. Other types of algorithms such as LARS jl8j are based on 



homotopy methods, and compute the whole solution path; i.e., they find the solution to (1.7) for 
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all values of the regularization parameter A and, in doing so, find the solution to the constrained 



problem (1.3). These methods do not scale well with problem size, however, especially when the 



solution is not that sparse. 

Turning to the Dantzig selector, solution algorithms are scarce. The standard way of solving 
(1.2) is via linear programming techniques 14 since it is well known that it can be recast as a 



linear program 15 . Typical modern solvers rely on interior-point methods which are somewhat 



problematic for large scale problems, since they do not scale well with size. Another way of solving 



(1.2) is via the new works 29 42 , which use homotopy methods inspired by LARS to compute the 



whole solution path of the Dantzig selector. These methods, however, are also unable to cope with 
large problems. Another alternative is adapting SPGL1 to this setting, but this comes with the 
caveat that it does not handle slight variations as discussed above. 



Finally, as far as the mixed norm problem (1.4) is concerned, we are not aware of efficient 



solution algorithms. One can always recast this problem as a second-order cone program (SOCP) 
which one could then solve via an interior-point method; but again, this is problematic for large- 
scale problems. 

1.3 Our approach 

In this paper, we develop a template for solving a variety of problems such as those encountered thus 
far. The template proceeds as follows: first, determine an equivalent conic formulation; second, 
determine its dual; third, apply smoothing; and fourth, solve using an optimal first- order method. 



1.3.1 Conic formulation 

In reality, our approach can be applied to general models expressed in the following canonical form: 



minimize 
subject to 



fix) 

A{x) + b e K. 



(1.9) 



The optimization variable is a vector x € M n , and the objective function / is convex, possibly 
extended-valued, and not necessarily smooth. The constraint is expressed in terms of a linear 
operator A : M. n —> M. m , a vector b £ M m , and a closed, convex cone /C C ]R m . We shall call a model 



of the form (1.9) that is equivalent to a given convex optimization model V a conic form for V . 

The conic constraint A(x)+b £ fC may seem specialized, but in fact any convex subset of W 1 may 
be represented in this fashion; and models involving complex variables, matrices, or other vector 
spaces can be handled by defining appropriate isomorphisms. Of course, some constraints are more 
readily transformed into conic form than others; included in this former group are linear equations, 
linear inequalities, and convex inequalities involving norms of affine forms. Thus virtually every 
convex compressed sensing model may be readily converted. Almost all models admit multiple 
conic forms, and each results in a different final algorithm. 



For example, the Dantzig selector (1.2) can be mapped to conic form as follows: 



\x\\ 1} A{x) ->■ (A* Ax, 0), b^(-A*y,5), K -> C[ 



(1.10) 



where £^ is the epigraph of the norm: = {(y,t) £ 



pn+1 . 



\y\\oo < t}. 
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1.3.2 Dualization 



The conic form (1.9) does not immediately lend itself to efficient solution using first-order methods 
for two reasons: first, because / may not be smooth; and second, because projection onto the set 
{x | A{x)-\-b G /C}, or even the determination of a single feasible point, can be expensive. We propose 
to resolve these issues by solving either the dual problem, or a carefully chosen approximation of 
it. Recall that the dual of our canonical form (1.9) is given by 



maximize 
subject to 



ff(A) 
A G /C*, 



(1.11) 



where g(\) is the Lagrange dual function 

5(A) ^ 

and K,* is the dual cone defined via 

K* = {A G 



inf C(x, A) = inf f(x) 

X X 



{X,A(x) + b), 



: (A, a;) > for all x G £}. 



The dual form has an immediate benefit that for the problems of interest, projections onto 
the dual cone are usually tractable and computationally very efficient. For example, consider the 
projection of a point onto the feasible set {x : \\Ax — y\\2 < e} of the LASSO, an operation which 
may be expensive. However, one can recast the constraint as A{x) + b £ JC with 



A(x) -)■ (Ax,0), 



/C 



■"2 ) 



(1.12) 



where is the second order cone £™ 
(Cf )* = C 



{(y,t) G R m+1 : \\y\\ 2 < t}. This cone is self dual, i.e., 
2 , and projection onto is trivial: indeed, it is given by 

\y,t), \\y\\ 2 < t, 

c{y, \\yh), -\\yh < * < Wvh, 
1(0,0), t<-\\ y \\ 2 , 



\vh + t 

%l|2 



(1.13) 



And so we see that by eliminating the affine mapping, the projection computation has been greatly 
simplified. Of course, not every cone projection admits as simple a solution as (1.13); but as we 
will show, all of the cones of interest to us do indeed. 



1.3.3 Smoothing 

Unfortunately, because of the nature of the problems under study, the dual function is usually 
not differentiable either, and direct solution via subgradient methods would converge too slowly. 
Our solution is inspired by the smoothing technique due to Nesterov [37]. We shall see that if one 
modifies the primal objective f(x) and instead solves 

minimize f^x) = f(x) + fid(x) ^ ^ 

subject to A(x) + b G /C, 

where d{x) is a strongly convex function to be defined later and /i a positive scalar, then the dual 
problem takes the form 

maximize g M (A) , . 

subject to A G JC*, 
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where is a smooth approximation of g. This approximate model can now be solved using first- 
order methods. As a general rule, higher values of [i improve the performance of the underlying 
solver, but at the expense of accuracy. Techniques such as continuation can be used to recover the 
accuracy lost, however, so the precise trade-off is not so simple. 

In many cases, the smoothed dual can be reduced to an unconstrained problem of the form 

maximize -j sm (z) — h(z), (1-16) 

with optimization variable z £ M m , where <7 sm is convex and smooth and h convex, nonsmooth, and 



possibly extended- valued. For instance, for the Dantzig selector (1.2), h(z) = 5\\z\\\. As we shall 
see, this so-called composite form can also be solved efficiently using optimal first-order methods. 
In fact, the reduction to composite form often simplifies some of the central computations in the 
algorithms. 

1.3.4 First-order methods 

Optimal first-order methods are proper descendants of the classic projected gradient algorithm. 



For the smoothed dual problem (1.15), a prototypical projected gradient algorithm begins with a 



point Aq £ /C*, and generates updates for k = 0, 1, 2, ... as follows: 



A fc+ i <- argmin ||A fc + tkVg^Xk) - A|| 2 , (1.17) 
Ae/C* 



given step sizes {t k }. The method has also been extended to composite problems like (1.16) 



38 , 46 49 ; the corresponding iteration is 



z k+ i <- argmin g sm {z k ) + {Vg sm (z k ),z - z k ) + ^-\\z - z k \\ 2 + h(z). (1-18) 
y 

Note the use of a general norm || • || and the inclusion of the nonsmooth term h. We call the 



minimization in (1.18) a generalized projection, because it reduces to a standard projection (1.17) if 
the norm is Euclidean and h is an indicator function. This generalized form allows us to construct 
efficient algorithms for a wider variety of models. 

For the problems under study, the step sizes {t k } above can be chosen so that e-optimality (that 



is, sup AgA; . 5^(A) — <?^(Afc) < e) can be achieved in 0(1/ e) iterations 36 . In 1983, Nesterov reduced 
this cost to 0(\/yfe) using a slightly more complex iteration 

Afc+i <r- argmin \\v k + t k Vg tl {v k ) - A|| 2 , v k+1 «- A fc+1 + a k {\ k+1 - X k ), (1.19) 
Ae/c* 

where vq = \q and the sequence {a k } is constructed according to a particular recurrence relation. 
Previous work by Nemirosvski and Yudin had established 0(1/ y/e) complexity as the best that 
can be achieved for this class of problems [33], so Nesterov's modification is indeed optimal. Many 
alternative first-order methods have since been developed [2j[30j|35j[37j|38j|46], including methods 
that support generalized projections. We examine these methods in more detail in ^5] 

We have not yet spoken about the complexity of computing or g sra and their gradients. For 
now, let us highlight the fact that V<7 At (A) = — „4(x(A)) — b, where 

x(X) = argmin C fM (x, A) = arg min f(x) + fid(x) — {A(x) + b, A), (1-20) 

X X 
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and d(x) is a selected proximity function. In the common case that d(x) 



X — Xq 



1 2 , the structure 



of (1.20) is identical to that of a generalized projection. Thus we see that the ability to efficiently 
minimize the sum of a linear term, a proximity function, and a nonsmooth function of interest is 
the fundamental computational primitive involved in our method. Equation (1.20) also reveals how 
to recover an approximate primal solution as A approaches its optimal value. 



1.4 Contributions 

The formulation of compressed sensing models in conic form is not widely known. Yet the convex 
optimization modeling framework CVX [23] converts all models into conic form; and the compressed 
sensing package £i-Magic |14| converts problems into second-order cone programs (SOCPs). Both 
systems utilize interior-point methods instead of first-order methods, however. As mentioned above, 



the smoothing step is inspired by 37 , and is similar in structure to traditional Lagrangian aug- 
mentation. As we also noted, first-order methods have been a subject of considerable research. 

Taken separately, then, none of the components in this approach is new. However their combi- 
nation and application to solve compressed sensing problems leads to effective algorithms that have 
not previously been considered. For instance, applying our methodology to the Dantzig selector 
gives a novel and efficient algorithm (in fact, it gives several novel algorithms, depending on which 
conic form is used). Numerical experiments presentedlater in the paper show that one can solve 
the Dantzig selector problem with a reasonable number of applications of A and its adjoint; the 
exact number depends upon the desired level of accuracy. In the case of the LASSO, our approach 
leads to novel algorithms which are competitive with state-of-the-art methods such as SPGL1. 

Aside from good empirical performance, we believe that the primary merit of our framework 
lies in its flexibility. To be sure, all the compressed sensing problems listed at the beginning of this 
paper, and of course many others, can be solved via this approach. These include total-variation 
norm problems, ^-analysis problems involving objectives of the form ||Wa;||i where W is neither 
orthogonal nor diagonal, and so on. In each case, our framework allows us to construct an effective 
algorithm, thus providing a computational solution to almost every problem arising in sparse signal 
or low-rank matrix recovery applications. 

Furthermore, in the course of our investigation, we have developed a number of additional tech- 
nical contributions. For example, we will show that certain models, including the Dantzig selector, 
exhibit an exact penalty property: the exact solution to the original problem is recovered even when 
some smoothing is applied. We have also developed a novel continuation scheme that allows us to 
employ more aggressive smoothing to improve solver performance while still recovering the exact 
solution to the unsmoothed problem. The flexibility of our template also provides opportunities to 
employ novel approaches for controlling the step size. 



1.5 Software 

This paper is accompanied with a software release [5], including a detailed user guide which gives 
many additional implementation details not discussed in this paper. Since most compressed sens- 
ing problems can be easily cast into standard conic form, our software provides a powerful and 
flexible computational tool for solving a large range of problems researchers might be interested in 
experimenting with. 

The software is not a single, monolithic solver; rather, it is a suite of programs and routines 
designed to serve as building blocks for constructing complete algorithms. Roughly speaking, we 
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can divide the routines into three levels. On the first level is a suite of routines that implement 
a variety of known first-order solvers, including the standard projected gradient algorithm and 
known optimal variants by Nesterov and others. On the second level are wrappers designed to 
accept problems in conic standard form (1.9) and apply the first-order solvers to the smoothed 
dual problem. Finally, the package includes a variety of routines to directly solve the specific 
models described in this paper and to reproduce our experiments. 

We have worked to ensure that each of the solvers is as easy to use as possible, by providing 
sensible defaults for line search, continuation, and other factors. At the same time, we have sought 
to give the user flexibility to insert their own choices for these components. We also want to provide 
the user with the opportunity to compare the performance of various first-order variants on their 
particular application. We do have some general views about which algorithms perform best for 
compressed sensing applications, however, and will share some of them in $6} 



1.6 Organization of the paper 

In §[2j we continue the discussion of conic formulations, including a derivation of the dual conic 
formulation and details about the smoothing operation. Section [3] instantiates this general frame- 
work to derive a new algorithm for the Dantzig selector problem. In £j4j we provide further selected 
instantiations of our framework including the LASSO, total-variation problems, ^-analysis prob- 
lems, and common nuclear-norm minimization problems. In £[5j we review a variety of first-order 
methods and suggest improvements. Section [6] presents numerical results illustrating both the em- 
pirical effectiveness and the flexibility of our approach. Section [7] provides a short introduction 
to the software release accompanying this paper. Finally, the appendix proves the exact penalty 
property for general linear programs, and thus for Dantzig selector and basis pursuit models, and 
describes a unique approach we use to generate test models so that their exact solution is known 
in advance. 



2 Conic formulations 



2.1 Alternate forms 



In the introduction, we presented our standard conic form (1.9) and a specific instance for Dantzig 



selector in (1.10). As we said then, conic forms are rarely unique; this is true even if one disregards 



simple scalings of the cone constraint. For instance, we may express the Dantzig selector constraint 
as an intersection of linear inequalities, —51 ^ A*(y— Ax) ^ 51, suggesting the following alternative: 



Fill) 



A(x) 



-A* A 
A* A 



51 
51 



A*y 
A*y 



/C 



a 2?i 



(2.1) 



We will return to this alternative later in £3.5 In many instances, a conic form may involve the 



manipulation of the objective function as well. For instance, if we first transform (1.2) to 



minimize t 
subject to ||x||i < t 

WA^y-Ax)^ < 5, 



then yet another conic form results: 

f(x,t)^t, A(x,t) ->■ (x,t, A*Ax,0), 



(2.2) 
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where C{ is the epigraph of the l\ norm, £" = {(y, t) G R n+1 : ||y||i < t}. 

Our experiments show that different conic formulations yield different levels of performance 
using the same numerical algorithms. Some are simpler to implement than others as well. Therefore, 
it is worthwhile to at least explore these alternatives to find the best choice for a given application. 



2.2 The dual 



To begin with, the conic Lagrangian associated with (1.9) is given by 

C(x,X) = f(x)-{X,A(x) + b), (2.3) 

where A £ M m is the Lagrange multiplier, constrained to lie in the dual cone tC* . The dual function 
g : R m -)• (MU -co) is, therefore, 

g(X) = inf C(x, A) = -f*(A*(X)) - (b, A). (2.4) 

X 

Here, A* : M m — > W 1 is the adjoint of the linear operator A and /*:!"-> (1U +co) is the convex 
conjugate of / defined by 

f*(z) = swp(z,x) - f(x). 



Thus the dual problem is given by 



maximize -f*(A*(X)) - (b,X) 
subject to A G IC* . 



(2.5) 



Given a feasible primal/dual pair (x, A), the duality gap is the difference between their respective 
objective values. The non-negativity of the duality gap is easily verified: 

f(x) - g(X) = f{x) + /* (A*(X)) + (b, A) > (x, A*(\)) + (b, A) = (A(x) + b,X)> 0. (2.6) 

The first inequality follows from the definition of conjugate functions, while the second follows 
from the definition of the dual cone. If both the primal and dual are strictly feasible — as is the 
case for all problems we are interested in here — then the minimum duality gap is exactly zero, and 
there exists an optimal pair (x*,A*) that achieves f(x*) = g(X*) = C(x*,X*). It is important to 
note that the optimal points are not necessarily unique; more about this in ^2.4. But any optimal 
primal/dual pair will satisfy optimality conditions 

A(x*) + b£tC, A*e/C\ (A(x*) + b, A*) = 0, A*(X*) e df(x*), (2.7) 

where df refers to the subgradient of /. 



2.3 The differentiable case 

The dual function is of course concave; and its derivative (when it exists) is given by 

Vg{X) = -A(x(X)) -b, x{X) e argmin£(x, A). (2.8) 

X 

It is possible that the minimizer x(A) is not unique, so in order to be differentiable, all such 
minimizers must yield the same value of — ^4(x(A)) — b. 
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If g is finite and differentiable on the entirety of JC* , then it becomes trivial to locate an 
initial dual point (e.g., A = 0); and for many genuinely useful cones K*, it becomes trivial to 



project an arbitrary point A G W 71 onto this feasible set. If the argmin calculation in (2.8) is 



computationally practical, we may entertain the construction of a projected gradient method for 



solving the dual problem (2.5) directly; i.e., without our proposed smoothing step. Once an optimal 



dual point A* is recovered, an optimal solution to the original problem (1.9) is recovered by solving 
x* G arguing C(x, A*). 

Further suppose that / is strongly convex; that is, it satisfies for some constant mf > 0, 

/((l - a)x + ax') < (1 - a)f{x) + af(x') - m f a(l - a)\\x - a/|||/2 (2.9) 

for all x,x' G dom(/) and < a < 1. Then assuming the problem is feasible, it admits a unique 
optimal solution. The Lagrangian minimizers x(X) are unique for all A G M n ; so g is differentiable 
everywhere. Furthermore, [37] proves that the gradient of g is Lipschitz continuous, obeying 

HVflr(A') - V<7(A)|| 2 < m/WAfWX' - A|| 2 , (2.10) 

where ||.4|| = sup|i 3 .|i 2=1 ||^4(x)||2 is the induced operator norm of A. So when / is strongly convex, 
then provably convergent, accelerated gradient methods in the Nesterov style are possible. 

2.4 Smoothing 

Unfortunately, it is more likely that g is not differentiable (or even finite) on all of fc* . So we 
consider a smoothing approach similar to that proposed in [37] to solve an approximation of our 
problem. Consider the following perturbation of ( |1.9[ ): 

minimize f^x) = f(x) + fj,d(x) _ _ 

subject to A(x) + b G K \ ■ > 

for some fixed smoothing parameter [i > and a strongly convex function d obeying 

d{x) > d{x a ) + \\\x - x Q \\ 2 (2.12) 



for some fixed point xq G W 1 . Such a function is usually called a proximity function. 

The new objective / M is strongly convex with mj = fx, so the full benefits described in { 2.3 
apply. The Lagrangian and dual functions becom 



now 



C^x, A) = f{x) + nd{x) - (A, A(x) + b) (2.13) 
g^(X) 4 inf C^x, X) = -(/ + »d)*(A*(X)) - (6, A). (2.14) 

One can verify that for the affine objective case f(x) = (cq,x) + do, the dual and smoothed dual 
function take the form 

g(X) = d -I {oy (A*(X)-co)-(b,X), 
g^X) = d - /j,d* (fi" 1 (A* (X) - c )) - (b, A), 



One can also observe the identity g M (A) = sup z — /*(^l*(A) — z) — fid*(n 1 z) — {b, A). 
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where /{g} is the indicator function of the set {0}; that is, 



J {0}(2/) 



0, y = 0, 
+oo, y ^ 0. 



The new optimality conditions are 

Aix^ + beJC, Xfj, G /C*, (^(x Ai ) + 6,A M ) = 0, X(A M ) -^Vd(^) e (2.15) 



Because the Lipschitz bound (2.10) holds, first-order methods may be employed to solve ( 2.11| ) 



with provable performance. The iteration counts for these methods are proportional to the square 
root of the Lipschitz constant, and therefore proportional to /z" 1 / 2 . There is a trade-off between 
the accuracy of the approximation and the performance of the algorithms that must be explored]^] 

For each [i > 0, the smoothed model obtains a single minimizer x^; and the trajectory traced 
by x M as (jl varies converges to an optimal solution x* = lim M _ > o + x^. Henceforth, when speaking 
about the (possibly non-unique) optimal solution x* to the original model, we will be referring to 
this uniquely determined value. Later we will show that for some models, including the Dantzig 
selector, the approximate model is exact: that is, x^ = x* for sufficiently small but nonzero \x. 

Roughly speaking, the smooth dual function is what we would obtain if the Nesterov smooth- 



ing method described in 37 were applied to the dual function g. It is worthwhile to explore how 



things would differ if the Nesterov approach were applied directly to the primal objective f(x). Sup- 
pose that f(x) = \\x\\i and d(x) = jll^lll- The Nesterov approach yields a smooth approximation 
fj} whose elements can be described by the formula 

[/?(*)].= sup zxi - = If I*;'-"' * = l,2,...,n. (2.16) 

|*|<i " - \x\ > fi, 

Readers may recognize this as the Huber penalty function with half-width fj,; a graphical comparison 
with fp is provided in Figure [TJ Its smoothness may seem to be an advantage over our choice 
fn(x) = \\x\\i + glMli) but the difficulty of projecting onto the set {x \ A(x) + b G K} remains; so 
we still prefer to solve the dual problem. Furthermore, the quadratic behavior of around xi = 
eliminates the tendency towards solutions with many zero values. In contrast, f^ix) maintains the 
sharp vertices from the i\ norm that are known to encourage sparse solutions. 

2.5 Composite forms 

In most of the cases under study, the dual variable can be partitioned as A = (2, r) G M™" m x W n 
such that the smoothed dual g^(z, r) is linear in r; that is, 

= -ffsm(z) - (vh,t) (2.17) 
for some smooth convex function <7 sm and a constant vector G M. m . An examination of the 



Lagrangian Cn (2.13) reveals a precise condition under which this occurs: when the linear operator 



A is of the form A{x) —> (A(x),0mxi), as seen in the conic constraints for the Dantzig selector 



(1.10) and LASSO (1.12). If we partition b = (b,b T ) accordingly, then evidently = b T . 



4 In fact, even when the original objective is strongly convex, further adding a strongly convex term may be 
worthwhile to improve performance. 
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Figure 1: The original objective f(x) — \x\ (blue), our modification (red), and Nesterov's smoothing 
(green) . 

Under such conditions, it is more natural to work with A, b, and b T directly, and exploit some 
useful simplifications. Specifically, let us define the function 

h : R m ~ ,fh — )> (IR U +oo), h(z) = inf{(6 r ,r) | (z,t) G JC*}. (2.18) 

Then the dual problem reduces to a nonsmooth unconstrained maximization 

maximize g^(z) = -g sm (z) - h(z). 

The gradient of g sm is Vg sm (z) = A(x(z)) + b, where x(z) is the minimizer of a reduced Lagrangian 

C^(x,z) = f(x) + nd(x) - (z,A(x) + b). (2.19) 

2.6 Projections 

A standard gradient projection step for the smoothed dual problem is 



Afc+i <- argmin||A- A fc - tfcV^(A fe )|| 2 . 
Ae/C* 



(2.20) 



For the composite version of the same problem, the corresponding generalized projection is 

z k +i <- argmin g sra (z k ) + {Vg sra (z k ), z - z k ) + ^-\\z - z k \\ 2 + h(z). (2.21) 

z K 

Integrating the definition of Vg sm (z) into ( 2.21[ ) and simplifying yields a two-sequence recursion: 

x k <- argmin/(x) + fid(x) - (A*(z k ),x) 



Zk+i <- argminh(z) + ^-\\z - z k \\ 2 + {A(x k ) + b, z). 



(2.22) 
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Note the similarity in computational structure of the two formulae. This similarity is even more 
evident in the common scenario where d(x) = h\\x — xo\\ 2 for some fixed xq £ M 71 . 

Let S be the matrix composed of the first m — fh rows of the m x m identity matrix, so that 



EA = z for all A = (z, r). Then (2.21) can also be written in terms of fC* and g^: 



z k+1 <r- E • argmax^(Afc) + (Vg^(\ k ), A - X k ) - ^r||S(A - A fc )|| , (2.23) 
where X k = (z k , h(z k )). If m = and the norm is Euclidean, then E = / and the standard projection 



(2.20) is recovered. So (2.21) is indeed a true generalization, as claimed in { 1.3.4 

The key feature of the composite approach, then, is the removal of the linear variables r from 
the proximity term. Given that they are linearly involved to begin with, this yields a more accurate 
approximation of the dual function, so we might expect a composite approach to yield improved 
performance. In fact, the theoretical predictions of the number of iterations required to achieve 
a certain level of accuracy are identical; and in our experience, any differences in practice seem 
minimal at best. The true advantage to the composite approach is that generalized projections 
more readily admit analytic solutions and are less expensive to compute. 

3 A Novel Algorithm for the Dantzig selector 

We now weave together the ideas of the last two sections to develop a novel algorithm for the 



Dantzig selector problem (1.2). 



3.1 The conic form 



We use the standard conic formulation (1.9) with the mapping (1.10) as discussed in the introduc- 
tion, which results in the model 

minimize ||x||i , . 

subject to (A*(y - Ax), 8) G ( ' 

where £™ is the epigraph of the norm. The dual variable A, therefore, will lie in the dual cone 

(3.2) 



(/2'JJ* = the epigraph of the l\ norm. Defining A = (z,t), the conic dual (2.5) is 

maximize —Ig 00 (—A*Az) — (A*y,z) — 5t 



subject to {z,t) G £f, 

where /* = Ii x is the indicator function of the H.^ norm ball as before. Clearly the optimal value 
of r must satisfy ||z||i = rj^jso eliminating it yields 

maximize —Ig ac (—A*Az) — (A*y,z) — 8\\z\\x. 



In both cases, the dual objectives are not smooth, so the smoothing approach discussed in {2.4 will 
indeed be necessary. 

5 We assume S > here; if 5 = 0, the form is slightly different. 
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3.2 Smooth approximation 

We augment the objective with a strongly convex term 

minimize IMIi + fJ>d(x) 
subject to (A* \y - Ax), 5) € K = 

The Lagrangian of this new model is 

£^(x; z,t) = \\x\\i + fid(x) - (z,A*(y - Ax)) - 5t. 

Letting x(z) be the unique minimizer of C^x; z, r), the dual function becomes 

g„(z,T) = ||x(z)||i +fid(x(z)) - (z,A*{y- Ax(z))) -t5. 



(3.3) 



Eliminating r per £2.5 yields a composite form gn(z) = —g S m( z ) — h(z) with 

g sm (z) = -||x(jk)||i - nd{x{z)) + (z,A*(y - Ax(z))), h{z) = S\\z\\i. 

The gradient of g sm is Vg sm (z) = A*(y- Ax(z)). 

The precise form of x(z) and V5f sm depend of course on our choice of proximity function d(x). 
For our problem, the simple convex quadratic 

d(x) = l\\x- x \\l, 

for a fixed center point xq £ W 1 , works well, and guarantees that the gradient is Lipschitz continuous 
with a constant at most /x -1 ||j4*^4|| 2 . With this choice, x(z) can be expressed in terms of the soft- 
thresholding operation which is a common fixture in algorithms for sparse recovery. For scalars x 
and r > 0, define 



SoftThreshold(x, t) = sgn(x) • max{|x| — r, 0} = < 



x + t, x < — r, 
0, |x| < r, 

X — T, X > T. 



When the first input a; is a vector, the soft-thresholding operation is to be applied componentwise. 
Armed with this definition, the formula for x(z) becomes 

x(z) = SoftThreshold(x - ^ 1 A*Az,^ 1 ). 

If we substitute x(z) into the formula for g sm _(z) and simplify carefully, we find that 

g sm (z) = SoftThresholdGuxo - A*Az, 1)||| + (A*y,z) + c, 

where c is a term that depends only on constants [i and xq. In other words, to within an additive 
constant, g S m( z ) is a smooth approximation of the nonsmooth term Ii oo (—A*Az) + {A*y,z) from 
(3.2), and indeed it converges to that function as [i — > 0. 
For the dual update, the generalized projection is 

Zk+i <- argmin5f sm (2; fc ) + (Vg sm (z k ), z - z k ) + ±-\\z - z k \\l + 5\\z\\i. (3.4) 

z k 

A solution to this minimization can also be expressed in terms of the soft thresholding operation: 

z k+1 <- SoftThreshold(z fe - t k A*(y - Ax(z k )),t k 5). 
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3.3 Implementation 



To solve the model presented in §3.2[ we considered first-order projected gradient solvers. After 
some experimentation, we concluded that the Auslender and Teboulle first-order variant |2,46| is a 
good choice for this model. We discuss this and other variants in more detail in so for now we 
will simply present the basic algorithm in Listing [T] below. Note that the dual update used differs 
slightly from (3.4) above: the gradient Vg sm is evaluated at y kl not z k , and the step size in the 
generalized projection is multiplied by O^ 1 . Each iteration requires two applications of both A and 
A* , and is computationally inexpensive when a fast matrix- vector multiply is available. 



Listing 1 Algorithm for the smoothed Dantzig selector 

Require: zq,xq £ W 1 , fi > 0, step sizes {t k } 
1: 6»o <- 1, v <- Zq 
2: for k = 0,1,2,... do 

3: y k 4- (1 - 9 k )v k + 9 k z k 

4: x k <— SoftThreshold(xo — n .A*Ay k , /U -1 )- 

5: z k+1 4- SoftThreshold(z fc - e~H k A*(y - Ax k ), B~H k 5) 

6: v k+1 4- (1 - 6 k )v k + 9 k Z k+ i 

7: 9 k+1 4- 2/(1 + (1+4/^)1/2) 
8: end for 



It is known that for a fixed step size t k = t < fi/\\A* A\\2, the above algorithm converges 
in the sense that g^{z*) — g^{z k ) = 0{k~ 2 ). Performance can be improved through the use of a 
backtracking line search for z k , as discussed in £ |5.3| (see also (3)). Further, fairly standard arguments 



in convex analysis show that the sequence {x k } converges to the unique solution to (3.3) 



3.4 Exact penalty 



Theorem 3.1 in |:8 can be adapted to show that as /x — > 0, the solution to (3.3) converges to a 



solution to (1.2). But an interesting feature of the Dantzig selector model in particular is that if 
fi < fj,Q for ^0 sufficiently small, the solutions to the Dantzig selector and to its perturbed variation 
(3.3) coincide; that is, x* = 2*. In fact, this phenomenon holds for any linear program (LP). 



Theorem 3.1 (Exact penalty) Consider an arbitrary LP with objective (c,x) and having an 
optimal solution (i.e., the optimal value is not —00) and let Q be a positive semidefinite matrix. 
Then there is a fiQ > such that ifO < \x < hq, any solution to the perturbed problem with objective 
(c,x) + \[i{x — xq,Q(x — xq)) is a solution to LP. Among all the solutions to LP, the solutions to 
the perturbed problem are those minimizing the quadratic penalty. In particular, in the (usual) case 
where the LP solution is unique, the solution to the perturbed problem is unique and they coincide. 

The theorem is proved in Appendix [A) There are also two related results in recent literature: 



a proof of the special case of noiseless basis pursuit [51], and a more general proof 21 that allows 



for a range of penalty functions. Our theorem is a special case of 21 , but we present the proof 
since it uses a different technique than |21| , is directly applicable to our method in this form, and 
provides useful intuition. The result, when combined with our continuation techniques in £5.5, is 
also a new look at known results about the finite termination of the proximal point algorithm when 
applied to linear programs [6l[40] . 
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Figure 2: Demonstration of the exact penalty property for the Dantzig selector. For /x < //o 
solution to the smoothed model coincides with the original model to solver precision. 



0.025, the 



As a consequence of the theorem, the Dantzig selector and noiseless basis pursuit, which are 
both linear programs, have the exact penalty property. To see why it holds for the Dantzig selector, 
recast it as the LP 

minimize 

subject to — u < x < u 

-51 < A*(y-Ax) < 81, 

with optimization variables (x,u) £ IR 2n . Then the perturbation ^n\\x — xo\\ 2 corresponds to a 
quadratic perturbation with a diagonal Q >z obeying Qa = 1 for 1 < i < n and Qu = for 
n + 1 < i < 2n. 

An illustration of this phenomenon is provided in Figure [2j For this example, a Dantzig selector 
model was constructed for a data set built from a DCT measurement matrix of size 512 x 4096. 
The exact solution x* was constructed to have 60 dB of dynamic range and 129 nonzeros, using 
the techniques of Appendix [Bj The figure plots the smoothing error \\x* — x^W as a function of fi; 
below approximately \x ~ 0.025, the error drops rapidly to solver precision. 

Unfortunately, compressed sensing models that are not equivalent to linear programs do not in 
general have the exact penalty property. For instance, it is possible to construct a counter-example 
for LASSO where the constraint is of the form \y — Ax\ < e. However, the result still suggests that 
for small e and ji, the LASSO and its smoothed version are similar. 



3.5 Alternative models 

As previously mentioned, different conic formulations result in different algorithms. To illustrate, 



consider the first alternative (2.1 ) proposed in ^2.1 which represents the Dantzig selector constraint 



via linear inequalities. The conic form is 

minimize 
subject to 



x \\i 

51 + A*{y 
51 - A*(y 



Ax) 
Ax) 



D 2n 



(3.5) 



16 



The dual variable A = (Ai, A2) must also lie in M+ n . The Lagrangian of the smoothed model is 

A) = ||x||i + lfi\\x - X0II2 - (Ai,<5l + A*{y - Ax)) - (A 2 , 61 - A*(y - Ax)) 

and its unique minimizer is given by the soft-thresholding operation 

x(A) = Soft Threshold^ - n~ 1 A*A(X 1 - X 2 ),^ 1 ). 

We cannot eliminate any variables by reducing to composite form, so we stay with the standard 
smoothed dual function ff^(A) = inf^ L^x; A), whose gradient is 



-51- A*(y- Ax(X)) 
-81 + A* \y - Ax{X)) 



V<7m(A) = 

The dual update is a true projection 

Afe+i = argmin-c^(A fc ) - (Vc^(A fc ), A - A fe ) + ^t^WX - A fe ||| 



= argmin||A - X k - t k Vg^(X k )\\2 
whose solution is simply the non-negative portion of a standard gradient step: 
A fc+ i = Pos(A fe + t k X7g^{X k )), [Pos(»]j 



(3.6) 



0, «j < 0. 



In order to better reveal the similarities between the two models, let us define z = Ai — A2 and 
f = l*(Ai + A2). Then we have ||z||i < f, and the Lagrangian and its minimizer become 

A) = ||x||i + \n\\x - x Q \\l - (z, A*(y - Ax)) - Sf, 
x{X) = SoftThreshold(x - fi~ 1 A*Az, A* -1 ); 

which are actually identical to the original norm-based model. The difference lies in the dual 
update. It is possible to show that the dual update for the original model is equivalent to 

= S-argmin-^(A fe )- (V^(A fc ),A-A fe ) + it fc 1 ||S(A-A fc )||2 (3.7) 

for E = [+7, — /]. In short, the dual function is linear in the directions of Ai + A2, so eliminating 
them from the proximity term would yield true numerical equivalence to the original model. 

4 Further Instantiations 

Now that we have seen the mechanism for instantiating a particular instance of a compressed 
sensing problem, let us show how the same approach can be applied to several other types of models. 
Instead of performing the full, separate derivation for each case, we first provide a template for our 
standard form. Then, for each specific model, we show the necessary modifications to the template 
to implement that particular case. 
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4.1 A generic algorithm 



A careful examination of our derivations for the Dantzig selector, as well as the developments in 
§2.6[ provide a clear path to generalizing Listing [T] above. We require implementations of the 
linear operators A and its adjoint A*, and values of the constants b, b T ; recall that these are the 
partitioned versions of A and b as described in {2.5 We also need to be able to perform the 
two-sequence recursion (2.22), modified to include the step size multiplier 9 k and an adjustable 
centerpoint xq in the proximity function. 

Armed with these computational primitives, we present in Listing [2] a generic equivalent of the 
algorithm employed for the Dantzig selector in Listing [TJ It is important to note that this particular 
variant of the optimal first-order methods may not be the best choice for every model; nevertheless 
each variant uses the same computational primitives. 



Listing 2 Generic algorithm for the conic standard form 
Require: z$,x$ G W 1 , /U > 0, step sizes {t k } 

1: 6»o <- 1, v <- Zq 

2: for k = 0,1,2,... do 

3: y k +- (1 - 9 k )v k + B k z k 

4: x k «- argmm x f(x) + (jd(x - x ) - (A*{y_ k ),x) 

5: z k+ i «- argmin z h(z) + ^\\z - z k \\ 2 + {A{x k ) + b, z) 

6: v k+1 k )v k + 9 k z k+ i 

7: 9 k+1 +- 2/(1 + (1 + 4/02)1/2) 

8: end for 



In the next several sections, we show how to construct first-order methods for a variety of 
models. We will do so by replacing lines 4-5 of Listing [2] with appropriate substitutions and 
simplified expressions for each. 

4.2 The LASSO 

The conic form for the smoothed LASSO is 

minimize [|a;||i + ^fi\\x — »o[|i f4 11 

subject to (y — Ax, e) G 

where is the epigraph of the Euclidean norm. The Lagrangian is 

C^x; z,t) = \\x\\i + \n\\x - xq\\\ — (z,y — Ax) - er. 

The dual variable A = (z,t) is constrained to lie in the dual cone, which is also Eliminating r 
(assuming e > 0) yields the composite dual form 

maximize inf x ||x||i + — xolli — ( z i V ~ Ax) — e||z||2- 

The primal projection with f(x) = \\x\\i is the same soft-thresholding operation used for the Dantzig 
selector. The dual projection involving h(z) = e||^||2, on the other hand, is 

z k+ i = argmine||2;||2 + ^\\z - z k \\\ + (x, z) = Shrink^ - 6 k 1 t k x,6~ 1 t k e), 
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where x = y — Ax k and Shrink is an ^-shrinkage operation 



Shrink^, t) = max{l — t/||z||2, 0} • z = 
The resulting algorithm excerpt is given in Listing [3j 



o, IMl2<t, 

{l-t/\\z\\ 2 )-z, ||z|| 2 >*. 



Listing 3 Algorithm excerpt for LASSO 



4: x k <— SoftThreshold(a;o — \i 1 A*y k ,fi 1 ) 
5: z k+1 <- Shrink^ - 0^Hk(y - Ax k ), 0^ x t k e) 



4.3 Nuclear-norm minimization 

Extending this approach to the nuclear-norm minimization problem 



(4.2) 



minimize ||AT||* 
subject to \\y — A{X)\\2 < e 

is straightforward. The composite smoothed dual form is 

maximize infx + fid(X) — (z,y — A(X)) — e\\z\\2, 

and the dual projection corresponds very directly to the LASSO. Choosing d{X) = ±\\X - X \\ 2 F 
leads to a primal projection given by the soft-thresholding of singular values: 

X k = SoftThresholdSingVal(X - H~ x JC(y h ),n~ x ). (4.3) 

The SoftThresholdSingVal operation obeys 

SoftThresholdSingVal(X, t) = U • SoftThreshold(E, t) • V* , 

where X = U"SV* is any singular value decomposition of Z, and SoftThreshold(S) applies soft- 
thresholding to the singular values (the diagonal elements of S). This results in the algorithm 
excerpt presented in Listing |4| 

Listing 4 Algorithm for nuclear-norm minimization (LASSO constraint) 

4: X k <- SoftThresholdSingVal(X - iT x A*{y k ), V' 1 ) 
5: z k+1 <- Shrink^ - 9^H k (y - A(X k )),e~H k e) 

Another constraint of interest is of Dantzig-selector type [IT] so that one is interested in 

minimize \\X\\* , , 

subject to \\A*(A(X) -y)\\< 5. ^ ' 

The cone of interest is, therefore, tC = {(X, t) : \\X\\ < t} and the dual cone is JC* = {(X, t) : 
\\X II* < t}. The derivation proceeds as before, and the composite dual problem becomes 

maximize inf x + \n\\X - X \\ 2 F - {Z,A*(y - A(X))) - 6\\Z\\*. 
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The gradient of the smooth portion has a Lipschitz continuous gradient with constant at most 
/x _1 ||.4*.4|| 2 , and singular value thresholding is now used to perform the dual projection. The 
resulting excerpt is given in Listing [5j 

Listing 5 Algorithm for nuclear-norm minimization (Dantzig-selector constraint) 

4: X k «- SoftThresholdSingVal(X - fx' 1 A* (A(Y k )) , ^) 

5: Z k+1 <- SoftThresholdSingVal(Z fc - &?t k A*(y - A(X k )),9^H k 5) 



4.4 ^x-analysis 

We are now interested in 



minimize ||Wx||i , 
subject to \\y — Ax\\2 < e, 



where the matrix W is arbitrary. This problem is frequently discussed in signal processing and 
is sometimes referred to as the method of l\- analysis. As explained in the introduction, this is a 
challenging problem as stated, because a generalized projection for f(x) = \\Wx\\i does not have 
an analytical form. 

Let us apply our techniques to an alternative conic formulation 

minimize t 

subject to ||Wx||i < t, 

\\y ~ Ax\\ 2 < e, 

where t is a new scalar variable. The dual variables are where 

lk (1) IU<rW, ||z( 2 )|| 2 < T ( 2 ), 

and the Lagrangian is given by 

C(x, t; Z W, r« z<?\ r^) = t - ( z W,Wx) - r«t - <z {2) , y - Ax) - er< 2 >. 

The Lagrangian is unbounded unless = 1; and we can eliminate in our standard fashion 
as well. These simplifications yield a dual problem 

maximize (y, z^) — e\\z^ ||2 
subject to A*zW _ W*z^ = 0, 

IN (1) ||oo<l- 

To apply smoothing to this problem, we use a standard proximity function d(x) = ~\\x — xq\\ 2 . 
(Setting = 1 causes t to be eliminated from the Lagrangian, so it need not appear in our 
proximity term.) The dual function becomes 

ffu(« (1) ,z (2) ) =inf ln\\x-x \\l - (z (1 \Wx) - (z^, y - Ax) - e\\z {2) II 

X 

and the minimizer x(z) is simply 

x(z) =x + l T 1 {W*z^ - A*z&). 
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Now onto the dual projection 

z k+ i = argmin e||z (2) || 2 + §t\\z - z k \\ 2 + (x, z), 

z: |k( 2 >||oo<l 

where x = (Wx(z), y — Ax{z)). This will certainly converge if the step sizes t k are chosen properly. 
However, if W and A have significantly different scaling, the performance of the algorithm may 

(i) 

suffer. Our idea is to apply different step sizes t k to each dual variable 



2 

z k+l = argmin e\\z^ || 2 + (x, z) + \B k J^)" 1 " 4° 111 

Z- lk (2) l|oc<l i= l 



in a fashion that preserves the convergence properties of the method. The minimization problem 
over z is separable, and the solution is given by 

= Trunc^ - ^Hfx^^H^) (4.6a) 
4 2) = Shrink(yf - ff^H^x^^H^e), (4.6b) 

where the truncation operator is given element-wise by 




Trunc(z,r) = sgn(z) • min{|z|, r} 

In our current tests, we fix = od, , where we choose a = ||W|| 2 /||.A|| 2 , or some estimate 
thereof. This is numerically equivalent to applying a single step size to a scaled version of the 
original problem, so convergence guarantees remain. In future work, however, we intend to develop 
a practical method for adapting each step size separately. 
The algorithm excerpt is given in Listing [6j 

Listing 6 Algorithm excerpt for ^i-analysis 

4: x k <-x + irHWyV>-A*yV t) ) 
K . 4'i <" TruncO/W - O^H^Wx^H^ 



z 



<- Shrink^ - 9^\y - Ax k ), O^e) 



4.5 Total-variation minimization 

We now wish to solve 



minimize |p||tv (a 7\ 

subject to \\y — Ax\\2 < e 



for some image array x 6 W 1 where ||x||tv is the total- variation introduced in 51.1 
actually cast this as a complex ^-analysis problem 

minimize ||Dx||i 
subject to \\y — Ax\\2 < e 



We can 
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where D : M™ 2 — > C^ n ^ 2 is a matrix representing the linear operation that places horizontal and 
vertical differences into the real and imaginary elements of the output, respectively: 

[Dxjij = (x i+ ij - xij) + \f^i ■ (scy+i - Xij), 1 < i < n, 1 < j < n. 

Writing it in this way allows us to adapt our ^i-analysis derivations directly. The smoothed dual 
function becomes 

gJzW,z^)=miU\\x-x \\ 2 2 - (z^,Dx) - (z™, y - Ax) - e\\z^ || 2 , 

X 

where z^ € W n is identical to the previous problem, and z^ G ci"-- 1 ) 2 satisfies H-z^Hoo < 1- 
Supporting a complex z^ requires two modifications. First, we must be careful to use the real- 
valued inner product 

(z^,Dx) 4 ^((z (1 YDx) = {^(D H z W )) T x 
Second, the projection requires a complex version of the truncation operation: 

[CTrunc(z,T)] k = min{l, r/\z k \} ■ z k = \ Zkl Zk \ ^ T ' 

\jZ k /\Z k \, \Z k \ > T. 

The algorithm excerpt is given in Listing [7} 

Listing 7 Algorithm excerpt for TV minimization 

4: x k <- xo + n~ l {^(D*y ( k ] ) ~ A*y ( k 2) ) 
5 . <- CTrunc(y« - rf^rf) 

• zfl <- Shrink(yf - fl- 1 ^ - Acjfe), ^ x 4 2) e) 



4.6 Combining £ x analysis and total-variation minimization 

We could multiply our examples indefinitely, and we close this section by explaining how one 



could solve the problem (1.4), namely that of finding the minimum of the weighted combination 



l + A| x[ tv subject to quadratic constraints. This problem can be recast as 



minimize 
subject to 



t + As 
||Wx||i < t 
\\Dx\\i < s 
\\Ax - y\\ 2 < e 



(4. 



and the strategy is exactly the same as before. The only difference with £4.4 and £4.5 is that the 



dual variable now belongs to a direct product of three cones instead of two. Otherwise, the strategy 
is the same, and the path is so clear that we prefer leaving the details to the reader, who may also 
want to consult the user guide which accompanies the software release [jjj. 
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5 Implementing first-order methods 



So far we have demonstrated how to express compressed sensing problems in a specific conic form 
that can be solved using optimal first-order methods. In this section, we discuss a number of 
practical matters that arise in the implementation of optimal first-order methods. This work 
applies to a wider class of models than those presented in this paper; therefore, we will set aside 
our conic formulations and present the first-order algorithms in their more native form. 

5.1 Introduction 

The problems of interest in this paper can be expressed in an unconstrained composite form 

minimize <fi(z) = g{z) + h(z), (5-1) 

where g, h : W 1 — > (MU +oo) are convex functions with g smooth and h nonsmooth. (To be precise, 
the dual functions in our models are concave, so we consider their convex negatives here.) Convex 
constraints are readily supported by including their corresponding indicator functions into h. 



First-order methods solve (5.1) with repeated calls to a generalized projection, such as 



z k+1 4- avgmmg(z k ) + (Vg(z k ),z- z k ) + ^-\\z - z k \\ 2 + h(z), (5.2) 

z " 

where || • || is a chosen norm and t k is the step size control. Proofs of global convergence depend 
upon the right-hand approximation satisfying an upper bound property 

4>{zk+i) < g(z k ) + {Vg(z k ),z k+ i - z k ) + 2^11^+1 - £fc|| 2 + h(z k+1 ). (5.3) 

This bound is certain to hold for sufficiently small t k ; but to ensure global convergence, t k must 
be bounded away from zero. This is typically accomplished by assuming that the gradient of g 
satisfies a generalized Lipschitz criterion, 

\\Vg(x) - Vg(y)\\* < L\\x - y\\ Vx, y G dom</>, (5.4) 



where || • ||* is the dual norm; that is, = sup{(z,w) \ \\z\\ < 1}. Then the bound ( |5.3[ ) is 

guaranteed to hold for any tk < L^ 1 . Under these conditions, convergence to e accuracy — that 



is, 4>(z k ) — inf z (j)(z) < e — is obtained in O(Lfe) iterations for a simple algorithm based on (5.2), 
and 0{^/Lfe) for the so-called optimal methods 0§[38|[46). These optimal methods vary the 



calculation (5.2) slightly, but the structure and complexity remain the same. 



5.2 The variants 

Optimal first-order methods have been a subject of much study in the last decade by many different 
authors. In 2008, Tseng presented a nearly unified treatment of the most commonly cited methods, 



and provided simplified proofs of global convergence and complexity 46 . The elegance of Tseng's 
effort seems underappreciated, possibly due to the fact that his premature passing delayed formal 
publication. His work has greatly eased our efforts to compare the performance of various algorithms 
applied to our conic models. 

We constructed implementations of five of the optimal first-order variants as well as a standard 
projected gradient algorithm. To simplify discussion, we have given each variant a 2-3 character 
label. Listing [8] depicts N07, a variation of the method described by Nesterov in (36|[38] . 
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Listing 8 Nesterov's 2007 algorithm (N07). 



Require: zq 6 dome/;, Lipschitz estimate L 
1: z <~ Z , 9 <-l 
2: for k = 0,1,2,... do 

3: y k «- (1 - 0jfc)2rjfc + 9 k z k 

4: z fc+1 <- argmin 2 (0£ £* =Q e^Vgfa), z) + §0£L||* - ^o|| 2 + /»(*) 
5: Zfc+i «- argmin 2 (V5f(y fc ),z) + |L||.z - y k \\ 2 + 
6: fc+1 <-2/(l + (l + 4/02)V2) 
7: end for 



The other variants can be described described simply by replacing lines 4-5 as follows. 



TS: Tseng's single-projection simplification of N07 |46|. 
4: z k+l ^aigm\n z {elY J k i=0 er l Vg{y i ),z) + \elL\\z-z4 2 + h{z) 

5: £fc+l (1 — k )z k + 9 k Z k+ \ 

LLM: Lan, Lu, and Monteiro's modification of N07 1 30 1 . 



4: z k+ i <- argmm z (V g(y k ),z) + \0 k L\\z - z k \\ 2 + h(z) 
5: z k+ i <- argmin^V g(y k ), z) + \L\\z - y k \\ 2 + h(z) 

AT: Auslender and Teboulle's method from [2]. 
4: z fc+ i <- argmin 2 (V5((y fc ),2;) + |# fc L||z - z k \\ 2 + 
5: <— (1 — k )z k + 9 k Z k+ \ 

N83: Nesterov's 1983 method (34f37j; see also FISTA [3]. 
4: z fc+ i <- argmin 2 (Vc/(y fc ),2) + - y fc || 2 + /i(z) 
5: Compute z k+ \ to satisfy z fc+ i = (1 - 9 k )z k + 9 k z k+1 . 



GRA: The classical projected gradient generalization. 

1 J 



4: z k+l <- argmm z (V g(y k ),z) + \L\\z - y k \\ 2 + h(z) 
5: z k +i «— z k+ i 

Following Tseng's lead, we have rearranged steps and renamed variables, compared to their original 
sources, so that the similarities are more apparent. This does mean that simpler expressions of 
some of these algorithms are possible, specifically for TS, AT, N83, and GRA. Note in particular 
that GRA does not use the parameter 9 k . 

Given their similar structure, it should not be surprising that these algorithms, except GRA, 
achieve nearly identical theoretical iteration performance. Indeed, it can be shown that if z* is an 



optimal point for (5.1), then for any of the optimal variants, 

<p(z k+1 ) - <j)(z*) < \L9 2 k \\z Q - z*\\ 2 < 2Lk- 2 \\z - z*\\ 2 . (5.5) 

Thus the number of iterations required to reach e optimality is at most [-y/2L/e||zo — z *\\ 2 ~\ (again, 
except GRA). Tighter bounds can be constructed in some cases but the differences remain small. 

Despite their obvious similarity, the algorithms do have some key differences worth noting. First 
of all, the sequence of points y k generated by the N83 method may sometimes lie outside of dom <j>. 
This is not an issue for our applications, but it might for those where g{z) may not be differentiable 
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everywhere. Secondly, N07 and LLM require two projections per iteration, while the others require 
only one. Two-projection methods would be preferred only if the added cost results in a comparable 
reduction in the number of iterations required. Theory does not support this trade-off, but the 



results in practice may differ; see {6.1 for a single comparison. 



5.3 Step size adaptation 

All of the algorithms involve the global Lipschitz constant L. Not only is this constant often difficult 
or impractical to obtain, the step sizes it produces are often too conservative, since the global 



Lipschitz bound (5.4) may not be tight in the neighborhood of the solution trajectory. Reducing L 



artificially can improve performance, but reducing it too much can cause the algorithms to diverge. 
Our experiments suggest that the transition between convergence and divergence is very sharp. 

A common solution to such issues is backtracking: replace the global constant L with a per- 
iteration estimate L\~ that is increased as local behavior demands it. Examining the convergence 



proofs of Tseng reveals that the following condition is sufficient to preserve convergence (see 46 
Propositions 1, 2, and 3): 



g{zk+i) < g{yk) + (Vg(yk),z k +i - yk) + \L k \ 



Zk+i ~ Vk\ 



(5.6) 



If we double the value of Lk every time a violation of (5.6) occurs, for instance, then Lk will 



satisfy Lj, > L after no more than \log 2 (L / Lq)~\ backtracks, after which the condition must hold 
for all subsequent iterations. Thus strict backtracking preserves global convergence. A simple 
improvement to this approach is to update L k with raecx.{2L k , L}, where L is the smallest value of 
Lk that would satisfy (5.6). To determine an initial estimate Lq, we can select any two points zq, z\ 
and use the formula 

Lo = l|Vff(zo) - Vg(zi)\\*/\\z - zi\\. 



Unfortunately, our experiments reveal that (5.6) suffers from severe cancellation errors when 
g(zk+i) ~ g{yk)i often preventing the algorithms from achieving high levels of accuracy. More 
traditional Armijo-style line search tests also suffer from this issue. We propose an alternative test 
that maintains fidelity at much higher levels of accuracy: 



\(yk ~ z k +i,V g{z k+ i) - Vg{yk))\ < lL k \\z k+ i - Vk 



(5.7) 



It is not difficult to show that (5.7) implies ( |5.6[ ), so provable convergence is maintained. It is a 
more conservative test, however, producing smaller step sizes. So for best performance we prefer 
a hybrid approach: for instance, use (5.6) when g(yk) — g(zk+\) > 75(^+1) for some small 7 > 0, 



and use (5.7) otherwise to avoid the cancellation error issues. 



A closer study suggests a further improvement. Because the error bound (5.5) is proportional 
to LkO^,, simple backtracking will cause it to rise unnecessarily. This anomaly can be rectified by 
modifying 9 k as well as during a backtracking step. Such an approach was adopted by Nesterov 
for N07 in |38| ; and with care it can be adapted to any of the variants. Convergence is preserved if 

L k+1 el +1 /{\ - e k+1 ) > L k el (5.8) 

(c.f. [46] , Proposition 1), which implies that the 6k update in Line 6 of Listing [8] should be 

6 k+l 4r- 2/(1 + (1 + AL k+1 /e 2 k L k ) 1/2 ). (5.9) 
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With this update the monotonicity of the error bound (5.5 ) is restored. For N07 and TS, the update 
for Zk+\ must also be modified as follows: 

z k+1 <- argmin^Lfc^to^O^V^y;),^} + \0 2 k L k \\z - z f + h(z). (5.10) 

z 

Finally, to improve performance we may consider decreasing the local Lipschitz estimate L k 
when conditions permit. We have chosen a simple approach: attempt a slight decrease of L k at 
each iteration; that is, L k = aL k _\ for some fixed a G (0, 1]. Of course, doing so will guarantee that 
occasional backtracks occur. With judicious choice of a, we can balance step size growth for limited 
amounts of backtracking, minimizing the total number of function evaluations or projections. We 
have found that a = 0.9 provides good performance in many applications. 

5.4 Linear operator structure 

Let us briefly reconsider the special structure of our compressed sensing models. In these problems, 
it is possible to express the smooth portion of our composite function in the form 

g(z)=g(A*(z)) + (b,z), 



where g remains smooth, A is a linear operator, and b is a constant vector (see {2.5 we have 
dropped some overbars here for convenience). Computing a value of g requires a single application 
of A*, and computing its gradient also requires an application of A. In many of our models, the 
functions g and h are quite simple to compute, so the linear operators account for the bulk of the 
computational costs. It is to our benefit, then, to utilize them as efficiently as possible. 



For the prototypical algorithms of Section 5.2, each iteration requires the computation of the 



value and gradient of g at the single point y k , so all variants require a single application each of A 
and A*. The situation changes when backtracking is used, however. Specifically, the backtracking 



test (5.6) requires the computation of g at a cost of one application of A*; and the alternative test 
(5.7) requires the gradient as well, at a cost of one application of A. Fortunately, with a careful 
rearrangement of computations we can eliminate both of these additional costs for single-projection 
methods TS, AT, N83, and GRA, and one of them for the two-projection methods N07 and LLM. 

Listing [9] depicts this more efficient use of linear operators, along with the step size adaptation 
approach, using the AT variant. The savings come from two different additions. First, we maintain 
additional sequences z A)k and z Ajk to allow us to compute y A ,k = A*{y k ) without a call to A*. 
Secondly, we take advantage of the fact that 

(Vk ~ z k+i, Vs>(z fc+ i) - Vg(y k )) = (yA,k ~ ZA,k+i,Vg(z Ayk+ i) - Vg(y A ,k)), 

which allows us to avoid having to compute the full gradient of g; instead we need to compute only 
the significantly less expensive gradient of g. 

5.5 Accelerated continuation 

Several recent algorithms, such as FPC and NESTA [4p7], have empirically found that continuation 
schemes greatly improve performance. The idea behind continuation is that we solve the problem of 
interest by solving a sequence of similar but easier problems, using the results of each subproblem 



to initialize or warm start the next one. Listing 10 below depicts a standard continuation loop 



26 



Listing 9 AT variant with improved backtracking. 



Require: z G dom0, L > 0, a € (0, 1], /3 e (0, 1) 
1: z «- Zo, ^o,^o <- ^*(^o), f-i = +oo, L_i = L 
2: for fc = 0,1,2,... do 
3: L fe <- aL fe _i 
4: loop 

5: 6 k <- 2/(1 + (1 + 4L jk /02_ 1 L fe _ 1 )V2) (0 Q 4 X) 

6: j/j. «- (1 - 6»fc)z fc + 6» fc Zfc, VA,k «- (1 - 0k)zA,k + k ZA,k 

7: 9fc <~ Vg(yj[,k), 9k <~ + & 

8: <- argmin^(5 fc ,z) + /i(z) + ^-||z - y fc || 2 , Z4 (fc +l <- „4*(z fc+1 ) 

9: Zfc+l «- (1 - ^fc) 2 ^ + O k z k +i, ZA,k+l <- (1 - 0k)zA,k + &kZA,k 

10: L <- 2|(y^ jfe - z^HbVs^fc+i) - ffjfc)| /||zjm-i ~ 2/fclli 

11: if L k > L then break endif 

12: Lfc <— m&x{L k / (3, L} 

13: end loop 

14: end for 



for solving the generic conic problem in (1.14) with a proximity function d(x) 



2 W x 



x 



We 



have used a capital X and loop count j to distinguish these iterates from the inner loop iterates 
x k generated by a first-order method. 



Listing 10 Standard continuation 

Require: Yq, fj,Q > 0, /? < 1 
1: for j = 0,1,2,... do 

2: Xj+i <- argmin4 (a . )+6e)C f{x) + ^\\x - Yj\% 
3: Y ■ . i <- .V ; . | or <- V) 

4: «- 
5: end for 



Note that Listing 10 allows for the updating of both the smoothing parameter //,■ and the 
proximity center Yj at each iteration. In many implementations, Yj = Xq and only fj,j is decreased, 
but updating Yj as well will almost always be beneficial. When Yj is updated in this manner, the 
algorithm is known as the proximal point method, which has been studied since at least [41]. Indeed, 
one of the accelerated variants [2] that is used in our solvers uses the proximal point framework to 
analyze gradient-mapping type updates. It turns out we can do much better by applying the same 
acceleration ideas already mentioned. 

Let us suggestively write 

h(Y)= mm f(x) + ^\\x-Y\\l (5.11) 

where \x > is fixed and C is a closed convex set. This is an infimal convolution, and h is known 
as the Moreau-Yosida regularization of / [28] . Define 

Xy = argmin/(x) + — \\x — y|||- (5-12) 



2 



The map Y i— > Xy is a proximity operator 32 . 
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We now state a very useful theorem. 



Theorem 5.1 The function h (5.11) is continuously differentiable with gradient 



V/i(Y) = n(Y-X Y ). 



(5.13) 



fx. Furthermore, minimizing h is equivalent 



The gradient is Lipschitz continuous with constant L 
to minimizing f(x) subject to x G C . 

The proof is not difficult and can be found in Proposition 1.2.2.4 and §XV.4.1 in [28j; see also 
exercises 2.13 and 2.14 in Yn, where h is referred to as the envelope function of /. The Lipschitz 
constant is fj, since Xy is a proximity operator P, and I— P is non-expansive for any proximity 
operator. 

The proximal point method can be analyzed in this framework. Minimizing h using gradient 
descent, with step size t = l/L = l/fi, gives 



Y, 



j'+i 



tVh(Yj) = Y 3 - 1 f i(Y j - X Yj ) = X Y] , 



which is exactly the proximal point algorithm. But since h has a Lipschitz gradient, we can use 
the accelerated first-order methods to achieve a convergence rate of 0(.T 2 ), versus 0(j 1 ) for 
standard continuation. In Listing 11 we offer such an approach using a fixed step size t = 1/L and 
the accelerated algorithm from Algorithm 2 in 46 . This accelerated version of the proximal point 



algorithm has been analyzed in 26 where it was shown to be stable with respect to inexact solves. 



Listing 11 Accelerated continuation 

Require: Yq, liq > 

1: X <- Y 

2: for j = 0, 1,2, . . . do 

3: X j+1 <- argmin^+^K f(x) + &-\\x - Yj\\ 2 2 
4: Y j+1 < .V,.| • (X j+1 - X, ) 
5: (optional) increase or decrease /ij 
6: end for 



Figure [3] provides an example of the potential improvement offered by accelerated continuation. 
In this case, we have constructed a LASSO model with a 80 x 200 i.i.d. Gaussian measurement 
matrix. The horizontal axis gives the number of continuation steps taken, and the vertical axis gives 
the error \\Xj — x*\\% between each continuation iterate and the optimal solution to the unsmoothed 
model. Both simple and accelerated continuation have been employed, using a fixed smoothing 
parameter fij = /xoj^] A clear advantage is demonstrated for accelerated continuation in this case. 

For models that exhibit the exact penalty property, updating the proximity center Yj at each 
iteration yields an interesting result. Examining the proof in Appendix |Aj we see that the property 
depends not on the size of \x but rather on the size of — Yj \\ , where is projection of Yj on the 
optimal set (e.g., if x* is unique, then x*j = x*). Therefore, for any positive /i, the exact penalty 
property will be obtained if Yj is sufficiently close to the optimal set. This has some obvious and 
very useful consequences. 

6 If jj, is fixed, this is not continuation per se, but we still use the term since it refers to an outer iteration. 
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»"* r Regular continuation (fixed u) 
: —Accelerated continuation (fixed u) 

o" 1 "^ — 1 1 1 — 

5 10 15 

outer iteration 

Figure 3: A comparison of simple continuation and accelerated continuation applied to a small-scale LASSO 
model. The horizontal axis gives the number of continuation steps taken; the vertical axis gives the error 
between the continuation solution and the original, unsmoothcd model. The error is normalized to 1 at 
iteration 1. 



The use of accelerated continuation does require some care, however, particularly in the dual 
conic framework. The issue is that we solve the dual problem, but continuation is with the primal 
variable. When Y is updated in a complicated fashion, as in Listing [TTJ or if fj, changes, we no 
longer have a good estimate for a starting dual variable Xj, so the subproblem will perhaps take 
many inner iterations. In contrast, if Y is updated in a simple manner and \x is fixed, as in Listing 



10 then the old value of the dual variable is a good initial guess for the new iteration, and this 
warm start results in fewer inner iterations. For this reason, it is sometimes advantageous to use 
Listing [10} 

Figure [4] shows an example where accelerated continuation does not offer a clear benefit. The 
model solved is the Dantzig selector, with a 64 x 256 partial DCT measurement matrix, and the 
exact solution (see Appendix [B]) has 50 nonzeros and a dynamic range of 68 dB which makes the test 
quite challenging for the minimum l\ solution is not the sparsest. Both fixed values of \x and both 
type of continuation are employed. The horizontal axis gives the total number of inner iterations, 
and the vertical axis gives the error \\xk — x*\\2 between the current iterate and the solution to the 
(unsmoothed) model. When no continuation is employed, the tradeoff between solution accuracy 
(small fj,) and the the number of iterations to converge (large [/,) is evident. Continuation helps 
drastically, and high levels of accuracy can be obtained quickly. In this example, both forms of 
continuation perform similarly. 

In our experience, accelerated continuation usually matches or outperforms regular continuation, 
but the matter requires further study. In particular, performance is affected by the stopping criteria 
of the inner iteration, and we plan to investigate this in future work. The experiments in Figure [4] 
decreased the tolerance by a factor of two every iteration, with the exception of the final iteration 
which used a stricter tolerance. 
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■no continuation, n/10 
-no continuation, \i 
-no continuation, 10 n 
regular continuation 
-accelerated continuation 




8000 

iterations 



Figure 4: Comparing fixed smoothing and continuation strategies for a Dantzig selector model. The 
horizontal axis gives the number of inner iterations of each first-order method. For the two continuation 
strategies, the circles depict the completion of each continuation step. 



5.6 Strong convexity 

Suppose that our composite function <j) is strongly convex; that is, 

0(#) > <P{ z o) + g^vIN ~~ 1 1 2 G dome/) 

for some fixed > 0. It is well known that a standard gradient method will achieve linear 
convergence in such a case. Unfortunately, without modification, none of the so-called optimal 
methods will do so. Thus in the presence of strong convexity, standard gradient methods can 
actually perform better than their so-called optimal counterparts. 

If the strong convexity parameter uia, < is known or can be bounded below, the N83 
algorithm can be modified to recover linear convergence (37 38 , and will recover its superior per- 



formance compared to standard gradient methods. The challenge, of course, is that this parameter 
rarely is known. In [38], Nesterov provides two approaches for dealing with the case of unknown 
mj,, one of which is readily adapted to all of the first-order variants here. That approach is a 
so-called restart method: the algorithm is restarted after a certain number of iterations, using the 
current iterate as the starting point for the restart; or equivalently, the acceleration parameter 9^ 
is reset to 0n = 1. In theory, the optimal number of iterations between restarts depends on m^/L^, 
but linear convergence can be recovered even for sub-optimal choices of this iteration count [25] . 

To illustrate the impact of strong convexity on performance, we constructed a strongly quadratic 
unconstrained function with = 0.07 and = 59.1, and minimized it using the gradient method 
(GRA), the AT first-order method with and without restart, and the N83 method tuned to the exact 
values of (m^L^). Backtracking was employed in all cases; we have verified experimentally that 
doing so consistently improves performance and does not compromise the exploitation of strong 
convexity. As can be seen, while AT without restart initially performs much better than GRA, 
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iterations 

Figure 5: First-order methods applied to a strongly convex problem: GRA, AT without restart, AT with 
restart, and N83 tuned with knowledge of the strong convexity parameter. 



the linear convergence obtained by GRA will eventually overtake it. The N83 method is clearly 
superior to either of these, achieving linear convergence with a much steeper slope. When restart is 
employed, AT recovers linear convergence, and achieves near parity with N83 when restarted every 
100 iterations (the optimal value predicted by 25 is every 112 iterations). 

The benefits of exploiting strong convexity seem clearly established. In fact, while proofs 
of linear convergence assume global strong convexity, our experiments in ^6] show that applying 
these methods to models with local strong convexity can improve performance as well, sometimes 
significantly. Unfortunately, the restart method requires manual tuning on a case-by-case basis 
to achieve the best performance. This is certainly acceptable in many applications, but further 
research is needed to develop approaches that are more automatic. 



6 Numerical experiments 

The templates we have introduced offer a flexible framework for solving many interesting but 
previously intractable problems; for example, to our knowledge, there are no first-order algorithms 
that can deal with complicated objectives like f(x) = ||Wx||i + ||a;||ry for a non-diagonal and 
non-orthogonal W. This section shows the templates in use to solve such real- world problems. It 
also describes some of the details behind the numerical experiments in previous sections. 



6.1 Dantzig selector: comparing first-order variants 



46 



there has been little focus on comparing the various accelerated 

Since 



Other than Tseng's paper 

methods. Tseng's paper itself presents few simulations that differentiate the algorithms 
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All variants, fixed and backtracking steps, no restart 



AT method, various restart strategies 
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calls to A and A* 
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Figure 6: Comparing first order methods applied to a smoothed Dantzig selector model. Left: comparing 
all variants using a hxed step size (dashed lines) and backtracking line search (solid lines). Right: comparing 
various restart strategies using the AT method. 



our software uses interchangeable solvers with otherwise identical setups, it is easy to compare the 
algorithms head-to-head applied to the same model. 

For this comparison, we constructed a smoothed Dantzig selector model similar to the one 
employed in §3.4| above. The model used a partial DCT measurement matrix of size 512 x 2048, a 
signal with 128 nonzero values, and an additive noise level of 30 dB SNR. The smoothing parameter 
was chosen to be [i = 0.25, and we then employed the techniques of Appendix [B] to perturb the 
model and obtain a known exact solution. This reference solution had 341 nonzeros, a minimum 
magnitude of 0.002 and a maximum amplitude 8.9. The smoothed model was then solved using 
the 6 first-order variants discussed here, using both a fixed step size of t = 1/L = fj,/\\A\\ 2 and our 
proposed backtracking strategy, as well as a variety of restart intervals. 

The results of our tests are summarized by two plots in Figure [6} The cost of the linear operator 
dominates, so the horizontal axes give the number of calls to either A or A* taken by the algorithm. 
The vertical axes give the relative error \\xk — x* \\/\\x*\\. Because this is a sparse recovery problem, 
we are also interested in determining when the algorithms find the correct support; that is, when 
they correctly identify the locations of the 341 nonzero entries. Therefore, the lines in each plot 
are thicker where the computed support is correct, and thinner when it is not. 

The left-hand plot compares all variants using both fixed step sizes and backtracking line search, 
but with no restart. Not surprisingly, the standard gradient method performs significantly worse 
than all of the optimal first-order methods. In the fixed step case, AT performs the best by a 
small margin; but the result is moot, as backtracking shows a significant performance advantage. 
For example, using the AT variant with a fixed step size requires more than 3000 calls to A or 
A* to reach an error of 10 ; with backtracking, it takes fewer than 2000. With backtracking, the 
algorithms exhibit very similar performance, with AT and TS exhibiting far less oscillation than 
the others. All of the methods except for GRA correctly identify the support (a difficult task due 
to the high dynamic range) within 1000 linear operations. 

The right-hand plot shows the performance of AT if we employ the restart method described in 
£5.6 for several choices of the restart interval. We observe significant improvements in performance, 



32 



revealing evidence of local strong convexity. A restart interval of 200 iterations yields the best re- 
sults; in that case, a relative error of 10 -4 is obtained after approximately 1000 linear operations, 
and the correct support after only a few hundred operations. The other variants (except GRA, 
which is unaffected by restart) show similar performance improvements when restart is applied, al- 
though the two-projection methods (N07 and LLM) take about 50% longer than the one- projection 
methods. 

Of course, care should be taken when applying these results to other contexts. For instance, 
the cost of the projections here is negligible; when they are are more costly (see, for instance, £6.4), 
two-projection methods (N07 and LLM) will be expected to fare worse. But even among Dantzig 
selector models, we found significant variations in performance, depending upon sparsity, noise 
level, and smoothing. For some models, the two-projection methods perform well; and in others, 
such as when we have local strong convexity, gradient descent performs well (when compared to 
the other algorithms without restart). Overall, it seems there is no best algorithm, but we choose 
the AT algorithm as our default since in our experience it is consistently one of the best and only 
requires one projection per iteration. 



6.2 LASSO: Comparison with SPGL1 

As mentioned in the introduction, there are numerous algorithms for solving the LASSO. Yet the 
algorithm produced by our dual conic approach is novel; and despite its apparent simplicity, it 
is competitive with the state of the art. To show this, we compared the AT first-order variant 
with SPGL1 [47] , chosen because recent and extensive tests in [4] suggest that it is one of the best 
available methods. 

The nature of the SPGL1 algorithm, which solves a sequence of related problems in a root- 
finding-scheme, is such that it is fastest when the noise parameter e is large, and slowest when e = 0. 
To compare performance in both regimes, we constructed two tests. The first is an "academic" 
test with an s-sparse signal and no noise; however, we choose s large enough so that the LASSO 
solution does not coincide with the sparse solution, since empirically this is more challenging for 
solvers. Specifically, the measurement matrix A is a 2 13 x 2 14 partial DCT, while the optimal value 
x* was constructed to have s = 2 12 nonzeros. The second test uses Haar wavelet coefficients from 
the "cameraman" test image (Figure [8] (a)) which decay roughly according to a power law, and 
adds noise with a signal-to- noise ratio of 30 dB. The measurement matrix is also a partial DCT, 
this time of size 0.3 • 2 16 x 2 16 . 

Figure [7] shows the results from both tests, each plot depicting relative error \\xk — x*||/||a;*|| 
versus the number of linear operations. We see that both methods achieve several digits of accuracy 
in just a few hundred applications of A and its adjoint. SPGL1 outperforms a regular AT solver 
in the left-hand "academic" test; however, AT with accelerated continuation solves the problem 
significantly faster than SPGL1. The noiseless case exploits our method's strength since the exact 
penalty property holds. 

For the wavelet test, SPGL1 outperforms our method, even when we use continuation. Although 
AT with continuation achieves high levels of accuracy in fewer than 1000 operations, other tests 
confirmed that SPGL1 is often a little better than our method, especially for large e. But the dual 
conic approach is competitive in many cases, and can be applied to a wider class of problem^} 

7 SPGL1 is also flexible in the choice of norm, but notably, it cannot solve the analysis problem due to difficulties 
in the primal projection. 
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Figure 7: Comparisons of the dual solver with SPGL1. The plot on the left involves a noiseless basis 
pursuit model, while the plot on the right represents a noisy image model. 



6.3 Wavelet analysis with total-variation 

The benefit of our approach is highlighted by the fact that we can solve complicated composite 
objective functions. Using the solver templates, it is easy to solve the £i-analysis and TV problem 
from £4.6 We consider here a denoising problem with full observations; i.e., A = I. Figure [8] (a) 
shows the original image Xq, to which noise is added to give an image y = xq + z with a signal-to- 
noise ratio of 20 dB (see subplot (b)). In the figure, error is measured in peak-signal-to- noise ratio 
(PSNR), which for an m x n 2 image x with pixel values between in [0, 1] is defined as 



PSNR(x) = 20 log 



10 



X — XoH-F 



where xq is the noiseless image. 

To denoise, we work with a 9/7 bi-orthogonal wavelet transform W (similar to that in JPEG- 
2000) with periodic boundary conditions (the periodicity is not ideal to achieve the lowest distor- 
tion) . A simple denoising approach is to hard-threshold the wavelet coefficients Wx and then invert 
with W~ 1 . Figure [8] (c) shows the result, where the hard-threshold parameter was determined ex- 
perimentally to give the best PSNR. We refer to this as "oracle thresholding" since we used the 



knowledge of xq to determine the threshold parameter. As is common with wavelet methods 44 
edges in the figure induce artifacts. 

Figures [8] (d) , (e) and (f) are produced using the solvers in this paper, solving 

minimize a||Wa;||i + j3\\x 
subject to \\Ax — y\\2 < e 



II j_ M 
IITV+ 2 



\x-y 



|2 
F 



(6.1) 



Figure[8](d) employs wavelet analysis only (a = = 0), (e) just TV [a = 0, /? = 1), and (f) both 
(a = 1, (3 = 5). For the best performance, the matrix A was re-scaled so that all the dual variables 
are of the same order; see |5| for further discussion of scaling. 



The Frobenius term in (6.1) is of course for smoothing purposes, and it is possible to minimize 



its effect by choosing fj, small or using the continuation techniques discussed. But for denoising, its 
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(a) Original (b) Noisy version (25.6 dB PSNR) 




(c) Wavelet thresholded (28.3 dB PSNR) (d) Wavelet I i-analysis (29.0 dB PSNR) 




(c) TV minimization (30.9 dB PSNR) (f) Wavelet analysis + TV minimization (31.0 dB PSNR) 

Figure 8: Dcnoising an n = 256 2 image. 
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presence makes little difference; in fact, it may give more visually pleasing results to use a relatively 
large /x. So to determine fi, we started with an estimate like 

H = max(a||Wy||i,/3||2/||TvO/c 

with c ~ 500 and then adjusted to give reasonable results. We ultimately employed /j = 1 for (d), 
fx = 50 for (e), and \i = 160 for (f). 

The wavelet analysis run took 26 iterations, and was complete in about 5 seconds. As shown in 
image (d), it produced boundary effects that are noticeable to the eye, and similar to those produced 
by thresholded image (c). The TV model (e) and TV with wavelets model (f) took 37 iterations (3 
seconds) and 30 iterations (8 seconds), respectively. Both produced better reconstructions, both by 
PSNR and by visual inspection. The additional wavelet analysis term in plot (f) offers only minimal 
improvement over TV alone, but this may be due to our simple choice of wavelet transform. For 
example, undecimated wavelets are common in denoising and may give better results, but our point 
here is simplicity and to point out the flexibility of the framework. 



6.4 Matrix completion: expensive projections 



We consider the nuclear-norm minimization problem (4.2 ) of a matrix X £ R n i xn 2 i n the conic dual 
smoothing approach. For matrix completion, the linear operator A is the subsampling operator 
revealing entries in some subset E C [n\\ x [n<2\. With equality constraints (e = 0) and Xq = 0, 
gradient ascent on the dual is equivalent to the SVT algorithm of |8|, a reference which also 



considered non-equality constraints, e.g., of the form (4.2). 

In addition to our own interest in this problem, one of the reasons we chose it for this article 
is that it differs from the others in one key respect: its computational cost is dominated by one of 
the projections, not by the linear operators. After all, the linear operator in this case is no more 
than a set of memory accesses, while the primal projection requires the computation of at least the 
largest singular values of a large matrix. As a result, the considerations we bring to the design of 
an efficient solver are unique in comparison to the other examples presented. 

There are a number of strategies we can employ to reduce the cost of this computation. The 
key is to exploit the fact that a nuclear-norm matrix completion model is primarily of interest when 
its optimal value X* is expected to have low rank [12] . Recall from § |4.3| that the update of X^ 
takes the form 

X k = SoftThresholdSingVal(A - lT 1 A*{\), lC 1 ). 

Our numerical experiments show that if ji is sufficiently small, the ranks r/% = rank(Vfc) will remain 
within the neighborhood of r* = rank(V*). In fact, with /i sufficiently small and Xq = 0, the rank 
grows monotonically. 

There are a variety of ways we can exploit the low-rank structure of the iterates X k . By storing 
Xk = C/fcSfcI4 in factored form, we can reduce the storage costs from 0{n\n2) to 0{rk{n\ + n2))- 
The quantity A(Xk), used in the computation of the gradient of the dual function, can be computed 



efficiently from this factored form as well. Finally, using a Lanczos method such as PROPACK 31 
the cost of computing the necessary singular values will be roughly proportional to r^. By combining 
these techniques, the overall result is that the cost of singular value thresholding is roughly linear 
in the rank of its result. These techniques are discussed in more detail in jjjj. 

Smaller values of /i, therefore, reduce the ranks of the iterates Xk, thereby reducing the com- 
putational cost of each iteration. This leads to a unique tradeoff, however: as we already know, 
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Figure 9: Noiseless matrix completion using various first-order methods. 



smaller values of \x increase the number of iterations required for convergence. In practice, we have 
found that it is indeed best to choose a small value of fi, and not to employ continuation. 

For our numerical example, we constructed a rank-10 matrix of size n\ = rii = 1000, and 
randomly selected 10% of the entries for measurement — about 5 times the number of degrees of 
freedom in the original matrix. We solved this problem using three variations of GRA and five 
variations of AT, varying the step size choices and the restart parameter. The smoothing parameter 
was chosen to be \i = 10 -4 , for which all of the methods yield a monotonic increase in the rank 
of the primal variable. The pure gradient method took about 1.8 minutes to reach 10 relative 
error, while the AT method without restart took twice the time; the error is in the Frobenius norm, 
comparing against the true low-rank matrix, which is the solution to the unperturbed problem. 

The results of our experiments are summarized in Figure [9} The horizontal axis now gives the 
number of SoftThresholdSingVal operations, a more accurate measure of the cost in this case; the 
cost of a SoftThresholdSingVal call is not fixed, but we observe that the rank of the iterates quickly 
reaches 10 and most SoftThresholdSingVal calls have roughly the same cost. The salient feature of 
this figure is that while AT initially outperforms GRA, the linear convergence exhibited by GRA 
allows it to overtake AT at moderate levels of precision. Employing a restart method with AT 
improves its performance significantly; and for a restart intervals of 50 iterations, its performance 
approaches that of GRA, but it does not overtake it. 

What accounts for this performance reversal? First, we observe that AT requires two projections 
per iteration while GRA requires only one. This difference is inconsequential in our other examples, 
but not for this one; and it is due to our use of the backtracking line search. Switching to a fixed 
step size eliminated the extra projection, but the overall performance suffered significantly. We 
hope to identify a new line search approach that avoids the added cost revealed here. 

Second, the similarities of Figures [5] and [9] suggest the presence of strong convexity. Using an 
estimate of the decay rate for gradient descent with step size t = 1/Lf, and comparing with known 
decay rate estimates [36] gives an estimate of rrif = 0.0024. For this value of mj (and with Lj = 1), 
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Figure 10: Noisy matrix completion on a non-low-rank matrix, using various first-order methods. 



the optimal restart number from 25 is about 80, which is consistent with the plot. It is easy 
to verify, however, that the smoothed dual function is not strongly convex. 

The results suggest, then, that local strong convexity is present. The authors in [25] argue 
that for compressed sensing problems whose measurement matrices satisfy the restricted isometry 
property, the primal objective is locally strongly convex when the primal variable is sufficiently 
sparse. The same reasoning may apply to the matrix completion problem; the operator .4**4 is 
nearly isometric (up to a scaling) when restricted to the set of low-rank matrices. The effect may 
be amplified by the fact that we have taken pains to ensure that the iterates remain low-rank. 
Overall, this effect is not yet well understood, but will also be explored in later work. 

We note that if the underlying matrix is not low-rank, then the local strong convexity effect is not 
present, and gradient descent does not outperform the AT method. In this case, our experiments 
also suggest that the restart method has little effect. Likewise, when inequality constraints are 
added we observe that the unusually good performance of gradient descent vanishes. For both the 
non-low-rank and noisy cases, then, it may be beneficial to employ an optimal first-order method 
and to use continuation. 



Figure 10 demonstrates these claims on a noisy matrix completion problem. We constructed a 
50 x 45 matrix of rank 20, sampled 67% of the entries, and added white noise to yield a 30 dB SNR. 
By using this small problem size, we are able to use CVX [23] to compute a reference solution; this 
took 5.4 minutes. The figure depicts three different approaches to solving the problem: GRA and 
AT each with fi = 5- 10 -4 and no continuation; and AT with \x = 10~ 2 and accelerated continuation. 
The restart approach no longer has a beneficial effect, so it is not shown. Each solver was run for 
500 iterations, taking between 2.5 and 4 seconds, and the plot depicts relative error versus the 
number of singular value decompositions. As we can see, the advantage of GRA has been lost, and 
continuation provides significant improvement: to achieve a relative error of 10 -2 , the continuation 
method required only 100 SVDs. 
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6.5 ^-analysis 



We present a brief example of solving the ^i-analysis problem (4.5 ), which is used in sparse recovery 
when the primal variable x is not sparse itself but rather sparse or compressible in some other 
domain. The domain may be well-known, such as frequency space (W is a DFT) or a wavelet 
basis (W is a wavelet transform), and for these cases W is invertible and even orthogonal and the 
problem may be solved using a method like SPGL1. However, these bases are often chosen out 
of convenience and not because they best represent the signal. More appropriate choices may be 
overcomplete dictionaries W G W xn with p 3> n, such as the undecimated wavelet transform, or 
the multilevel Gabor dictionary. The experiment below uses the Gabor dictionary with p = 28n. 
To solve the LASSO problem using a dictionary W, the two common approaches are analysis: 

minimize ||Wa;||i , ^\ 

subject to \\y — Ax\\2 < e 

and synthesis: 

minimize ||a||i , , 

subject to \\y - AW*a\\2 < e, ^ ' ' 

with decision variable a. Similarly, the Dantzig selector approach would have the same objectives 
but constraints of the form ||^4*(y — Ae)||oo < 5 (analysis) and ||^4*(y — .AW* a) ||oo < 5 (synthesis). 
When W is not orthogonal, the two approaches are generally different in non-trivial ways, and 



furthermore, the solutions to synthesis may be overly sensitive to the data 19 . 

The differences between analysis and synthesis are not very well understood at the moment for 
two reasons. The first is that the analysis problem has not been studied theoretically. An exception 
is the very recent paper [9] which provides the first results for ^-analysis. The second is that there 
are no existing efficient first-order algorithms to solve the analysis problem, with the exception of 
the recent NESTA (2] and C-SALSA [I] algorithms, which both work on the LASSO version, and 
only when AA* = I. 

With the dual conic approach, it is now possible to solve the smoothed analysis problem, and by 
using continuation techniques, the effect of the smoothing is negligible. To illustrate, we constructed 
a realistic test of the recovery of two radio-frequency radar pulses that overlap in time, as depicted 



in Figure 11 The first pulse is large, while the second pulse has 60 dB smaller amplitude. Both 
have carrier frequencies and phases that are chosen uniformly at random, and noise is added so 
that the small pulse has a signal-to- noise ratio of 0.1 dB. 

The signal is recovered at Nyquist rate resolution for 2.5 GHz bandwidth, and the time period 
is a little over 1600 ns, so that n = 8,192. The sensing matrix A is modeled as a block-diagonal 
matrix with ±1 entries on the blocks, representing a system that randomly mixes and integrates 
the input signal, taking 8 measurements every 100 ns, which is 12.5 x below the Nyquist rate. Thus 
A is a 648 x 8,192 matrix and W is a 228,864 x 8,192 matrix. We applied both the Dantzig selector 
and LASSO models to this problem. 

To solve the problems, we employed the AT variant with accelerated continuation. At each 
iteration, the stopping tolerance is decreased by a factor of 1.5, and the continuation loop is 
ended when it no longer produces significantly different answers, which is usually between 2 and 8 
iterations. The value of fi is set to 

n , \\Wxls\\2 (a a\ 



2 



\ X LS\ 
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Figure 11: Recovery of a small pulse 60 dB below a large pulse. The first plot employs a Dantzig selector, 
the bottom employs LASSO. 
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Table 1: Details on the simulation used in Figure 11 For both the Dantzig selector and LASSO versions, 



the algorithm was run with continuation and reweighting. The "cont." column is the number of continuation 
steps, the "iter." column is the number of iterations (over all the continuation steps), the "time" column is 
in seconds, and the "error" column is the relative 1-2, error. Each row is one solve of the Dantzig selector or 
the LASSO. 



where xls is the least-squares solution to Ax = y, which is easy to calculate since m n (and 
independent of p) . 

To enhance the results, we employed an outer reweighting loop [17], in which we replace the 
||Wx||i term with ||i?Wx||i for some diagonal weight matrix R. Each reweighting involves a full 
solve of either the Dantzig selector or the LASSO. This loop is run until convergence, which is 
typically about 2 to 5 iterations. The results are plotted in the frequency domain in Figure 11 The 
large pulse is the dominant spike at about 1.1 GHz, and it is easily recovered to very high precision 
(the relative £2 error is less than 5 • 10~ 3 ). The small pulse is at about 2.2 GHz, and because it has 
an SNR of 0.1 dB and the measurements are undersampled by a factor of 12, it is not possible to 
recover it exactly, but we can still detect its presence and accurately estimate its carrier frequency. 

Table [T] reports the computational results of the test, first solving either the Dantzig selector or 
the LASSO, and then taking a single reweighting step and re-solving. Each call of the algorithm 
takes about 100 iterations, and the Dantzig selector or the LASSO is solved in about 1 minute, 
which is impressive since, due to the extremely large size of W, this problem is intractable using 
an interior-point method. 



7 Software: TFOCS 

The work described in this paper has been incorporated into a software package, Templates for 
First-Order Conic Solvers (TFOCS, pronounced tee-fox), which will be made publicly available at 
|http : //tf ocs . Stanford. edu[ As its name implies, this package is a set of templates, or building 
blocks, that can be used to construct efficient, customized solvers for a variety of models. 

To illustrate the usage of the software, let us show how to construct a simple solver for the 
smoothed Dantzig selector model described in ^3] We will begin by assuming that we are given the 
problem data A, b, delta, and a fixed smoothing parameter mu. The basic solver templates require 
two functions to complete their work. The first function computes the value and gradient of g sm , 
the smooth component of the composite dual: 

function [ val, grad, x ] = g_dantzig( A, y, xO, mu, z ) 
x = SoftThreshold( xO - (1/mu) * A' * ( A * z ), 1/mu ); 
grad = A ' * ( y - A * x ) ; 

val = z' * grad - norm( x, l)-0.5*mu* norm( x - xO ) . " 2; 

The second function computes the generalized projection associated with the nonsmooth component 
of the dual, which is in this case h(z) = 5\\z\\\. 
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function [ z, h ] = h_dantzig( delta, z_old, grad, L ) 
z = SoftThreshold( z_old - (1/L)*grad, delta/L ); 
h = delta * norm( z, 1 ); 

Both of these functions depend on the soft-thresholding operator: 

function y = Sof tThreshold( x, t ) 

y = sign( x ) . * max( abs( x ) - t, ); 

Armed with these functions, the following code solves the smoothed Dantzig selector, using the 
Auslender/Teboulle first-order variant and the default choices for line search and stopping criteria. 

function x = Dantzig_smoothed( A, y, delta, mu ) 
[m,n] = size (A) ; 

xO = zeros(n,l); zO = zeros(n,l); 
g_sm = @(z) g_dantzig( A, y, xO, mu, z ); 
h = @(z,g,L) h_dantzig( delta, z, g, L ); 
[z,x] = solver_AT( g_sm, h, zO ); 

The last line of code calls the AT solver, which upon completion delivers its solution to both the 
smoothed dual and the primal problem. 

This simple solver is likely to be unnecessarily inefficient if A exhibits any sort of fast opera- 
tor structure, but this can be remedied rather simply. Note that the solver itself needs to have 
no knowledge of A or y above; its interaction with these quantities comes only through calls to 
g_dantzig. Therefore, we are free to rewrite g_dantzig in a more numerically efficient manner: 
for instance, if A is derived from a Fourier transforms, we may substitute fast Fourier transform 
operations for matrix-vector multiplications. This simple change will reduce the cost of each itera- 
tion from 0(mn) to O(nlogn), and the storage requirements from 0(n 2 ) to 0{n). For large-scale 
problems these savings can be quite significant. 

A further improvement in performance is possible through careful management of the linear 



operator calculations as described in £5.4 The TFOCS solver templates can perform this manage- 
ment automatically. To take advantage of it, we replace g_dantzig with two functions: one which 
implements the underlying linear operator, and one which implements the remainder of the smooth 
function. The details of how to accomplish this are best left to the user guide [5] . 

Of course, as part of the complete package, we have supplied a solver for the Dantzig selector 
that exploits each of these efficiencies, and others. Similar solvers have been created for the LASSO, 
TV, and other models discussed here. But the flexibility of the lower-level templates will allow these 
same efficiencies to be achieved for many models that we have not discussed here. In fact, evidently 



the templates are not restricted to our specific conic form ( 1.9 ) or its dual, and we hope the software 
finds application outside of the compressed sensing domain as well. 

With the software release is a detailed user guide [S] that covers the usage, and documents the 
most popular problem formulations and provide some examples. We refer the reader to the user 
guide for further software details. 

8 Conclusion 

We have developed a convenient framework for constructing first-order methods, which is flexible 
and handles a variety of convex cone problems, including problems that did not have efficient 
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algorithms before. On the implementation side, we have introduced ideas which lead to novel, 
stable, and efficient algorithms. When comparing our implementation on specific problems such as 
the LASSO, which have been intensively studied, our techniques appear surprisingly competitive 
with the state of the art. 

The templates from this paper are flexible in a manner that has not yet been seen in sparse 
recovery software. Not only are the solvers interchangeable and hence easy to benchmark one 
algorithm against another, but they work with a wide variety of formulations which we hope will 
greatly enable other researches. It is also our goal that the software will be easy to use for non- 
experts, and to this end our future research will be to improve the usability of the software and 
to provide sensible default parameters. One major topic to pursue is choosing the smoothing 
parameter fi. Further efforts will also be made to improve the line search to take advantage of 
strong convexity, to better manage expensive projections, and to use scaled norms so that dual 
variables are all on the same scale (this is an issue only for objectives with several additive terms, 
as in the TV with analysis problem in £ |4.6| and £6.3). 

A subsequent paper will cover these issues, as well as further investigating local strong convexity 
and how to take advantage of this in an optimal manner, and improving the accelerated continuation 
scheme by taking into account the inexact solves. The software and user guide [5] will be kept up- 
to-date and supported. 



A Exact Penalty 

The general approach proposed in this paper is to add a strongly convex term to the primal objective 
function in order to smooth the dual objective. We now draw a parallel with augmented Lagrangian 
and penalty function methods that eliminate constraints by incorporating a penalty term p into 
the objective, e.g., mm x: Ax=b f(x) becomes min x f(x) + Xp(Ax — b). For reasonable choices of p, 
the two problems become equivalent as A — > oo. Remarkably, for some special choices of p, the two 
problems are equivalent for a large but finite value of A. Usually, p is a non-differentiable function, 
such as p = || • ||i or p = \\ ■ W2 (not p = \\ ■ |||); see Bertsekas' book [7j for a discussion. 

Our approach uses a different type of perturbation, since our goal is not to eliminate constraints 
but rather to smooth the objective. But we use the term "exact penalty" because, for some 
problems, we have a similar result: the smoothed and unsmoothed problems are (nearly) equivalent 
for some /j, > 0. The result below also departs from traditional exact penalty results because our 
perturbation ^\\x — xq\\\ ^ s smooth. 



Below we present the proof of Theorem 3.1 which exploits the polyhedral constraint set of 
linear programs. 



Proof of Theorem 3.1 We consider the linear program 

minimize (c, x) 



(LP) 

subject to x G V, 



where V is a convex polyhedron, and its perturbed version 

minimize (c, x) H 
subject to x G V, 



minimize (c, x) + ^(x - x ,x - xq)q 
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where (x, y) = (x, Qy) for some positive semidefinite Q. Let E* be the solution set of (LP) (E* is a 
vertex or a face of the feasible polyhedron) and let x* be any point in E* such that (x — xq, x — xq)q 
is minimum. (Note that when Q is the identity, x* is the usual projection onto the convex set 
E*.) With ffi(x) = (c, x) + \p{x — xq, x — xo)q, it is of course sufficient to show that x* is a local 
minimizer of /„ to prove the theorem. Put differently, let x E V and consider x* + t(x — x*) where 
< t < 1. Then it suffices to show that lim t _ >0 + /^(x* + t(x — x*)) > f^(x*). We now compute 

/ M (x* + i(x - x*)) = (c, x* + t(x - x*)) + \p{x* + t(x - x*) - x , x* + t(x - x*) - x ) Q 

= Ufa*) + t(c, X - X*) + flt(x - X*, X* - X )q + 7j/Ut 2 (x - x*, x - x*)q. 

Therefore, it suffices to establish that for all x E V, 

(c, x — x*) + fj,{x — X*, X* — Xq)q > 0, 

provided p is sufficiently small. Now x £ V can be expressed as a convex combination of its extreme 
points (vertices) plus a nonnegative combination of its extreme directions (in case V is unbounded). 
Thus, if {vi} and {dj} are the finite families of extreme points and directions, then 

x = ^2 ^ + 5^ pjdj> 

» j 

with Aj > 0, Yli^i = 1) l°i — 0- The extreme directions — if they exist — obey (c,dj) > as 
otherwise, the optimal value of our LP would be — oo. Let / denote those vertices which are not 
solutions to (LP ), i.e., such that (c, vi — x) > for any x 6 E* . Likewise, let J be the set of extreme 
directions obeying (c, dj) > 0. With this, we decompose x — x* as 

x-x* = (y^ X i( v i ~ x *) + p i d 3) + (2 Xi ( Vi ~ X *^ + Y1 Pi d i) ■ 
i£i j£J iei jeJ 

It is not hard to see that this decomposition is of the form 

x - x* = a(x - x*) + ^i( v i ~ x *) + ^2 Pjdjj , 

iei jeJ 

where x £ E* and a is a nonnegative scalar. This gives 

(c,x - x*) = ^Xii^Vi - x*) + y^ j p j (c,d j ). 
iei jeJ 

Let oti > (resp. /3j > 0) be the cosine of the angle between c and V{ — x* (resp. c and dj). Then 
for i E I and j € J, we have on > and /3j > 0. We can write 

(c, x - x*) = ^2 ^i°H\\ c h\\vi -x*h + y^, Pjl3j\\ch\\dj\\2, 
iei jeJ 

which is strictly positive. Furthermore, 

(x - x*,x* - x ) Q = a(x- x*,x* - x ) Q + (^A^u; - x*) + s ^p j dj,x k - x ^) 

iei jeJ Q 



iei jeJ 

> -\\Q(x* - x )h(j2 x *W Vi ~ x *h + J2pj\\ d jh)- 

iei jeJ 
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The second inequality holds because x* minimizes (x — xq,x — xo)q over E*, which implies {x — 
x*,x* — xq)q > 0. The third is a consequence of the Cauchy-Schwartz inequality. In conclusion, 
with hq = fj,\\Q(x* — xo)||2) we have 

(c, x - x*) + n{x - x*,x* - x )q > }Xai\\c\\ 2 - fiQ)Xi\\vi - x*\\ 2 + } )(Pj\\c\\2 - Ho)Pj\\dj\\2, 

and it is clear that since mim, > and mim f3j > 0, selecting \x small enough guarantees that the 
right-hand side is nonnegative. 



B Creating a Synthetic Test Problem 

It is desirable to use test problems that have a precisely known solution. In some cases, such as 
compressed sensing problems in the absence of noise, the solution may be known, but in general this 
is not true. A common practice is to solve problems with an interior-point method (IPM) solver, 
since IPM software is mature and accurate. However, IPMs do not scale well with the problem 
size, and cannot take advantage of fast algorithms to compute matrix-vector products. Another 
disadvantage is that the output of an IPM is in the interior of the feasible set, which in most cases 
means that it is not sparse. 

Below, we outline procedures that show how to generate problems with known exact solutions 
(to machine precision) for several common problems. The numerical experiments earlier in this 



paper used this method to generate the test problems. This is inspired by 38 , but we use a 



variation that gives much more control over the properties of the problem and the solution. 



Basis pursuit. Consider the basis pursuit problem and its dual 

minimize ||^||i maximize (y, A) 

subject to Ax = y, subject to ||A*A||oo < 1- 

At the optimal primal and dual solutions x* and A*, the KKT conditions hold: 

Ax" = y p*Aloo < 1 (A*X*) T = sign(x^), 

where T = supp(x*). 

To generate the exact solution, the first step is to choose A and y. For example, after choosing 
A, y may be chosen as y = Ax for some x that has interesting properties (e.g., x is s-sparse or is 
the wavelet coefficient sequence of an image). Then any primal dual solver is run to high accuracy 
to generate x* and A*. These solutions are usually accurate, but not quite accurate to machine 
precision. 

The idea is that x* and A* are exact solutions to a slightly perturbed problem. Define T = 
supp(x); in sparse recovery problems, or for any linear programming problem, we have \T\ < m 
where m is the length of the data vector y. The matrix A is modified slightly by defining A <— AD 
where D is a diagonal matrix. D is calculated to ensure that || (DA* X*)^ ||oo < 1 and (DA*X*)t = 
sign(a;*). If the original primal dual solver was run to high accuracy, D is very close to the identity. 
In practice, we observe that the diagonal entries of D are usually within .01 of 1. 
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The primal variable x* is cleaned by solving Aj-x*^ = y; this is unique, assuming At has full 
column rank. If the original problem was solved to high accuracy (in which case D will have all 
positive entries) , then the cleaning- up procedure does not affect the sign of x* . The vectors x* and 
A* are now optimal solutions to the basis pursuit problem using A and y. 



The LASSO. A similar procedure is carried out for the LASSO and its dual given by 

minimize ||x||i maximize (y, A) — e||A||2 

subject to \\Ax — 2/ 1 1 2 < e, subject to || A*A||oo < 1. 

Let z be the variable such that y = Ax + z, then strong duality holds if 

((Ax*, A*) - + ((z, X*) - e||A|| 2 ) = 0. 

The operator A is chosen as before, so (Ax*, A*) — ||x*||i = 0. Thus z needs to satisfy (z, A*) — e|| A 1 1 2 , 
i.e., z = ey/\\y\\2- This means that x* and A* are optimal solutions to the LASSO with A and 
y = Ax* + z. 



Other problems. For basis pursuit and the LASSO, it is possible to obtain exact solutions to 
the smoothed problem (with d(x) = ^\\ x ~ ^olll)- For the Dantzig selector, an exact solution for 
the smoothed problem can also be obtained in a similar fashion. To find an exact solution to the 
unsmoothed problem, we take advantage of the exact penalty property from £3.4 and simply find 
the smoothed solution for a sequence of problems to get a good estimate of xo and then solve for 
a very small value of fi. 
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